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ABSTRACT 


From  theoretical  studies  of  two-phase  flow  of  gas -particle  mixtures  in 
rocket  nozzles,  the  governing  equations  for  two-dimensional  and  axisymmetrical 
flows  have  been  derived.  Characteristic  equations  have  been  put  into  finite  dif¬ 
ference  form  and  a  mathematical  model  has  been  developed  for  a  general  interior 
point  in  the  supersonic  flow  field.  The  subsonic  and  transonic  flow  fields  have 
been  investigated.  To  integrate  the  subsonic,  transonic,  and  supersonic  solutions, 
a  mathematical  model  will  be  programed  for  the  IBM  7090  computer  that  will  de¬ 
termine  nozzle  performance  for  given  nozzle  contour  and  motor  conditions. 
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I.  INTRODUCTION 

Liquid  or  solid  particles,  or  a  combination  of  the  two,  exist  in  the  exhaust 
products  of  solid  propellants  when  powdered  metals  are  used  as  fuel  additives. 
These  additives  are  included  in  the  propellant  formula  for  several  reasons,  the 
foremost  being  to  increase  the  energy  release  of  the  propellant  and  to  prevent 
combustion  instability. 

However,  several  important  deleterious  effects  occur  when  condensable 
products  are  formed  in  the  exhaust  gases  as  a  result  of  having  metal  additivies 
in  the  propellant.  Significant  performance  losses  occur  because  the  drag  of  the 
gas  on  the  particles  is  irreversible:  the  thermal-energy  content  of  the  particles 
transfers  to  the  gas  by  convection  and,  to  some  degree,  by  radiation,  which  re¬ 
quires  a  finite  amount  of  time.  Since  the  particles  pass  through  the  nozzle  very 
rapidly,  this  energy  transfer  does  not  have  time  to  take  place  completely,  and 
a  significant  amount  of  the  thermal  energy  of  the  particles  is  lost.  Impingement 
of  particles  on  the  nozzle  wall,  affecting  the  heat  transfer  rates  and  causing 
erosion  of  the  wall,  alters  the  internal  configuration  of  the  nozzle  and  results  in 
a  change  in  performance. 

To  minimize  these  effects,  it  is  necessary  to  predict  the  effects  of  various 
geometric  and  gas -dynamic  parameters  on  the  performance  of  a  gas -particle 
mixture.  The  objective  of  this  study  was  to  develop  an  analytical  technique  for 
determining  the  effects  of  particles  on  nozzle  performance  and  design. 

n.  SUMMARY 


This  report  presents  a  theoretical  discussion  of  the  two-phase  flow  of 
gas -particle  mixtures  in  rocket  nozzles  and  shows  the  derivation  of  the  govern¬ 
ing  equations  for  two-dimensional  and  axisymmetric  flows.  Characteristic 
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II,  Summary  (cont. ) 

equations  are  put  into  finite-difference  form.  A  numerical  solution  procedure, 
developed  for  a  general  interior  point  in  the  supersonic  flow  field,  has  been 
programed  for  the  IBM  7090  computer.  Correctness  of  the  computer  program 
has  been  verified  through  hand  calculation  of  the  numerical  solution  derived  for 
the  IBM  7090. 

A  method  has  been  developed  for  determining  the  subsonic  flow  field  on  the 
basis  of  a  one -dimensional  sink-flow  approximation.  The  transonic  flow  field  has 
been  approximated  by  a  power -series  solution  of  the  linearized  system  equations 
for  constant  kinetic  and  thermal  lags  of  the  particles . 

A  method  is  then  suggested  for  integrating  the  subsonic,  transonic,  and 
supersonic  solutions  into  a  computer  program  with  which  nozzle  performance  for 
given  chamber  conditions  and  nozzle  contour  can  be  determined. 

This  report  indicates  that  status  of  gas -particle -flow  studies  as  of  February 
1962.  Since  that  date  significant  progress  has  been  made  on  a  program  that  in¬ 
tegrates  the  subsonic,  transonic,  and  supersonic  solutions  into  one  IBM  7090  com¬ 
puter  program  with  which  overall  nozzle  performance  can  be  evaluated.  The  pro¬ 
gram  enables  study  of  several  particle  sizes  and  chemical  species.  The  thermody¬ 
namic  and  flow  properties  of  the  gas  and  particles  are  evaluated  as  a  function  of 
temperature.  Frozen  or  shifting  chemical  equilibrium  can  be  studied  by  specify¬ 
ing  the  appropriate  speed  of  sound  as  a  function  of  gas  temperature.  Studies  can 
bjs  made  at  any  point  in  the  nozzle  if  a  supersonic  starting  line  can  be  determined 
by  the  use  of  other  techniques.  Currently,  this  program  has  been  limited  to  con¬ 
ventional  de  Laval  nozzles. 

m.  TECHNICAL  DISCUSSION 

A.  EQUATIONS  GOVERNING  GAS-PARTICLE  MIXTURES 
1.  Basic  Assumptions 
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III,  A,  Equations  Governing  Gas  ^Particle  Mixtures  (corit. ) 

Derivation  of  the  equations  governing  gas -particle  mixtures 
is  given  in  Appendix  A.  These  equations  were  derived  on  the  basis  of  the  follow 
ing  assumptions: 


a.  The  particles  are  spherical  and  all  the  same  size,. 

b.  The  total  mass  of  the  gas -particle  mixture  is  constant. 

c.  The  total  energy  of  the  gas -particle  mixture  is  constant. 

d.  The  internal  temperature  of  the  particle  is  uniform,  and 
the  particle  specific  heat  is  constant. 

e.  The  gas  and  particles  exchange  thermal  energy  by  con¬ 
vection  only. 

f.  The  gas  obeys  the  perfecUgas  law,  has  a  constant 
molecular  weight,  and  has  constant  specific  heats. 

g.  All  external  forces  except  pressure  of  the  gas  and  drag 
of  the  particles  are  neglected. 

h.  The  gas  is  inviscid  except  for  the  drag  it  exerts  on  the 
particles. 

i.  The  particles  do  not  interact  with  each  other. 


j.  The  volume  occupied  by  the  particles  is  negligible. 
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III,  A,  Equations  Governing  Gas -Particle  Mixtures  (cont.  ) 


2.  Governing  Equations 


For  two-dimensional  (0  =  ^>0  ^  =0)  and  axisymmetric  (a  =  y>^y~ 
steady  flows,  the  governing  equations  are: 


p  (u  )  +  P  (v  )  +  u  (p  )  +  V  {p  ) 

g  X  ^g'  g'y  g  ^g  X  g  ^g'y 


y  o 


p  (u  )  +  p  (v  )  +  u  (p  ;  +  V  (p  )  ~~a 

p'x  ’^p'  p'y  p  p  X  p  ^P  y  y  a 


0  V 

LJe._E. 


p  r  u  (U  )  +  V  (u  )  1  +  Ap  (u  -u  )  +  {P„), 

^g  L  8  8  *  8  g  y  J.  P  g  P  g  > 


=  0 


p  Pulv)  +v{v)"l  +  Ap  (v  -V  )  +  (P  )  =0 

^g  L  8  g  X  g'  gVJ  V  g  g  y 

u  (P_)  +  V  (P  )  r  a^r  u  (p  )  +  '^o^Pa^vl  “ 
g  gx  g'  g'y  L  8  gx  g  gyj  I 


u  (u  )  +  V  (u  )  »  A(u  -  u  ) 

p'  p'x  p'  p'y  g  p' 


u  (v  )  +  V  (v  )  =  A(v  -  V  ) 

P  P  X  p'  p'y  '  g  P 


u  (h  )  +  V  (h  )  =-  4  AC(T  -  T  ) 

p'  p'x  p'  p'y  T  '  p  g 


B  =  0 


P  =  p  RT 
g  '^g  g 


P  g 


(1) 

(2} 

(3) 

(4) 

{5) 

(6) 

(7) 

(8) 

(9) 
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III,  Technical  Discussion  (cont.  ) 


B.  CHARACTERISTIC  EQUATIONS 


By  the  use  of  the  method  of  characteristics  as  discussed  in  Appendix 
B  and  developed  in  detail  in  Appendix  C,  the  characteristic  equations  for  the  above 
system  of  quasi-linear  partial  differential  equations  (equation  1  through  8)  were 
found  to  be: 


1 .  Along  Gas  Streamlines 


g 

Pg[Ug%  t  Vgdvjt  dP^  =  -  V  ^  ■  'p' 


dP 

g 


dp. 


8  Pd 


Ap  B  dx 


8  8 


{17} 

(18) 

(19) 


2.  Along  Gas  Mach  Lines 


2  2 
u  -a 
8 

(Ugdy.Vgdx)[APpBdx  -  UgdPjf  ’“pHy  -  (v^.Vp)dx]dx 


(20) 


(21) 
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III,  Technical  Discussion  (cont.  ) 

B.  CHARACTERISTIC  EQUATIONS 

By  the  use  of  the  method  of  characteristics  as  discussed  in  Appendix 
B  and  developed  in  detail  in  Appendix  C,  the  characteristic  equations  for  the  above 
system  of  quasi-linear  partial  differential  equations  (equation  1  through  8)  were 
found  to  be: 


1 .  Along  Gas  Streamlines 


(17; 

(18} 

(19) 


(Ugdy  -Vgdx)  [APpBdx  -  u^dP  J+  {ap^  [(u.^  -Up)dy  -  (v^  -Vp)dx]  dx 

+Pg[vgdu^-Ugdv^-ay-^(u^dy-Vgdx)]dx  +  dP^  dyj.*0  (21) 
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III,  B,  Characteristic  Equations  (cont.  ) 


3.  Along  Particle  Streamlines 


,  V 

dy  _  p 

dx  ~  u 
P 

u  du  =  A(u  -  u  )  dx 
P  P  g  P 


V  dv  =  A(v  -V  )  dy 
P  P  g  P 

u  dh  =  -  4  AC  (T  -  T  )  dx 
P  P  3  p  g 


( )  =(T  p  u 
P  y  ^P  P 


(^  )  =  -<yp  V 
p  X  ^p  p 


(22) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 


This  system  of  characteristics  is  totally  hyperbolic  if  the  flow 
is  supersonic  (M>1),  and  partially  hyperbolic  when  the  flow  is  subsonic.  Hence, 
the  flow  in  the  supersonic  portion  of  the  nozzle  may  be  obtained  by  use  of  a  nu¬ 
merical  procedure  based  on  the  above  system  of  characteristic  equations.  How¬ 
ever,  the  subsonic  and  transonic  flow  must  be  obtained  differently,  since  the 
characteristics  are  only  partially  hyperbolic  in  such  regions. 

C.  SUBSONIC  SOLUTION 

Because  of  the  elliptic  nature  of  the  governing  equations  for  gas- 
particle  mixtures  in  subsonic  flow,  a  solution  could  not  be  obtained  for  the  two- 
dimensional  and  axisymmetric  flow  equations.  In  approximation,  the  subsonic 
flow  was  assumed  to  be  one  dimensional.  The  following  system  of  equations  was 
then  derived  (Appendix  O),  subject  to  the  same  assumptions  as  for  the  two- 
dimensional  and  axisymmetric  cases. 
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RT 

g 


Tabulated, 


area  a  f(x) 


(35) 

(36) 

(37) 

(38) 
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III,  C,  Subsonic  Solution  (cont.  ) 

The  first  three  differential  expressions  can  be  solved  by  techniques 
such  as  the  Adams  or  Runge-Kutta  methods  once  starting  values  are  assumed  at 
some  location  where  the  flow  velocities  are  small  and  gas -particle  equilibrium 
may  be  assumed. 

Since  lines  of  constant  properties  in  the  subsonic  portion  of  a  de  Laval 
nozzle  are  concave  downstream,  the  area  in  these  equations  should  be  based  on  the 
spherical  sector  determined  by  the  distance  to  the  vertex  of  an  equivalent  cone 
which  is  tangential  to  the  wall  of  the  nozzle  at  each  point.  This  is  equivalent  to 
a  sink-flow  solution,  where  the  location  of  the  sink  varies  with  the  position  in 
the  nozzle  and  the  slope  of  the  nozzle  wall,  as  illustrated  in  Figure  1. 


Figure  1.  One -Dimensional  Sink-Flow  Model 

Thus,  the  subsonic  georhetry  would  be  specified  by  giving  the  wall  location  y  and 
slope  G  as  a.  tabulated  function  of  x.  The  area  as  a  function  of  x  is  determined  as 
follows; 
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1II>  C,  Subsonic  Solution  (cont.  ) 

X  =  X  -  R  (  1  -  cos  0 )  (40) 

area  =  ZirR^  (  1  -  cos  0  )  (41) 


By  using  the  above  relationship  between  x  and  area,  the  subsonic 
solution  may  be  carried  out  to  a  point  where  the  solution  begins  to  diverge 

because  of  the  ( — ^ — )  term.  To  perform  this  solution,  the  total  mass  flow  must 
M  -1 

be  estimated  and  must  subsequently  be  corrected  after  the  starting  line  in  the 
supersonic,  region  is  determined. 

w  =  w  +  w  (42) 

g  P 


The  gas  mass-flow  rate  may  be  estimated  from  the  chamber  stagnation  condi¬ 
tions  by  assuming  that  the  gas -particle  mixture  is  in  equilibrium  and  by  using 
the  constant-lag  one -dimensional  analysis  presented  in  Appendix  E  to  determine 
the  velocity  and  density  at  some  downstream  cross  section,  such  as  the  nozzle 
throat,  of  an  equivalent  perfect-gas  flow. 

Wg  =  PgUg  area  {43V 


(44) 


(45) 
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III,  C,  Subsonic  Solution  (cont.  ) 


2 

u 

g 


(area)  _  J_r  2 
(area)*  ^  ^ 


(V  +t) 


(46) 


(47) 


E  =  1  +  {k[v  (1-k)  +  k]+  (V  -1) 


(50) 


1 


If  the  nozzle  throat  is  chosen  as  the  reference  cross  sectioni  then  Ki  si 
and  (area)  ^  »  (area)  .  Equations  (43)  through  (46)  then  reduce  to 


Report  No.  0162-01TN-16 


III,  C,  Subsonic  Solution  (cont.  ) 


To  evaluate  the  above  expressions,  a  value  for  the  kinetic  lag,  K,  must  be  as¬ 
sumed.  With  some  experience,  a  reasonable  average  value  of  K  between  the 
entrance  and  throat  may  be  chosen.  However,  as  a  first  approximation, 


w 


K  may  be  chosen  as  unity.  Since  ^  is  known,  w  can  be  calculated  when 


is  known.  Hence,  all  the  parameters  necessary  for  the  subsonic  solution  can 
be  determined. 


D.  TRANSONIC  SOLUTION 


The  one -dimensional  sink  flow  solution  does  not  converge  near  Mach  1, 

and  the  method  of  characteristics  cannot  be  used  where  the  Mach  number  is  less 

than  1.  Thus,  a  different  approach  was  necessary  in  the  transonic  flow  regime. 

The  particle  kinetic  and  thermal  lags  were  assumed  to  be  constant,  and  the 

analysis  shown  in  Appendix  E  for  a  constant-lag  one -dimensional  gas -particle 

mixture  was  used  to  obtain  a  modified  specific  heat  ratio,  Y_.  and  Mach  number, 

S 

of  the  gas,  which  includes  the  effects  of  the  particles  on  the  gas.  The  Sauer 
transonic  flow  approximation  discussed  in  Appendix  F  was  then  used  to  obtain  the 
gas  velocity  components  in  the  supersonic  regime  where  the  Mach  number  is  only 
slightly  greater  than  1.  The  other  gas  and  particle  properties  can  then  be  found 
along  this  starting  line  by  the  use  of  the  following  procedure.  A  starting  line  hav¬ 
ing  all  properties  known  is  available  for  use  in  the  supersonic  characteristic 
solution. 


The  transonic  flow  regime  was  broken  into  two  regions  of  constant 
lag:  the  first  from  the  end  of  the  sink-flow  regime  to  the  throat,  and  the  second 
from  the  throat  to  the  starting  line.  This  is  illustrated  in  Figure  2. 
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III,  D,  Transonic  Solution  (cont.  ) 


This  model  permits  different  wall  radii  of  curvature  to  be  used  upstream, 
Pj,  and  downstream,  of  the  throat.  In  Appendix  F,  the  gas  flow  field 
near  Mach  1  is  described  by  the  following  relationships. 


(Vg  + 

u(x,y)-ax  + 

(Y  +  1)  ax 

*  flit,)  y* 


2 

-y 

(y+  1)^°^ 
zO+'b)"  i^~b) 


3 


y 


e  =  - 


2(3+6) 


W(yg+  1)  (i  +  6)y, 

V  o 


(1+6) 

“V(v  +  P.Vs 


(56) 

(57' 


(59) 
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III,  D,  Transonic  Solution  (cont.  ) 


where  6=0  for  two-dimensional  flow  and 6=  1  for  axisymmetric  flow,  and  u  and  v 
are  nondimensional  perturbation  velocities.  The  actual  velocity  components  of 
the  modified  gas  flow  are  give:-,  by  the  follov/ing  expressions  containing  u.  and  v. 


u  (x ,  y )  =  a  *  ( 1  +  u ) 

(60) 

V  (x,  y)  =  a*  V 

(61) 

- S—  R  T 

(Vg+  1) 

(62) 

The  actual  velocity  components  of  the  gas  in  the  gas -particle  mixture  are  de¬ 
termined  from  the  following  relationship,  wherein  E  is  defined  as  a  one - 


dimensional  flow  parameter  with  constant  lag; 
2  _ 


u  =  E  u 

i63) 

V  =  E  V 

g 

(64) 

The  value  of  K  (one-dimensional  particle  velocity  lag)  used  in  Region  1 
(Figure  2)  is  determined  from  the  sink -flow  solution  along  the  sink  flow  line. 
When  the  entire  gas -particle  flow  field  is  finally  determined  at  the  throat,  the 
value  of  K  at  the  throat  may  be  found.  An  average  value  on  K  between  the  sink 
flow  line  and  the  throat  can  then  be  established,  thus  determining  the  Region  1 
flow  field  more  accurately.  This  procedure  should  be  repeated  until  the  average 
value  of  K  no  longer  varies  from  trial  to  trial.  Region  2  can  be  solved  in  the 
same  manner  by  first  choosing  K  as  the  value  at  the  throat  from  the  Region  1 
solution.  When  the  starting  line  is  determined,  an  average  value  of  K  in 
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III,  D,  Transonic  Solution  (cont. ) 

Region  2  can  be  repeated  in  the  same  manner  as  for  Region  1.  When  the  flow 
fields  in  Regions  1  and  2  are  completely  determined,  the  gas  velocity  components 
on  the  supersonic  starting  line  are  determined. 

Next,  it  is  necessary  to  make  the  total  mass  flow  assumed  in  the  sink-flow 
solution  compatible  with  the  total  mass  flow  in  the  supersonic  region.  This  may 
be  done  by  integrating  the  total  mass  flow  across  the  starting  line  to  obtain  the 
supersonic  mass  flow.  The  assumed  total  mass  flow  in  the  subsonic  sink-flow 
solution  should  then  be  set  equal  to  the  calculated  total  mass  flow  in  the  super¬ 
sonic  region,  and  the  entire  subsonic  and  transonic  solutions  should  be  repeated 
to  obtain  a  new  starting  line  and  supersonic  total  mass  flow.  This  procedure 
should  be  repeated  until  the  subsonic  and  supersonic  total  mass  flows  become 
equal;  the  supersonic  solution  may  then  be  initiated. 

The  particle  velocity  components  in  the  transonic  region  may  be  found  by 
using  the  method  of  characteristics,  since  the  particle  streamlines,  which  are 
also  characteristics,  exist  in  subsonic  as  well  as  supersonic  flow.  By  choosing 
several  points  on  the  sink  flow  line  as  starting  points,  the  particle  properties 
along  the  streamlines  through  these  points  may  be  found  by  using  numerical 
techniques.  The  geometry  for  this  procedure  is  illustrated  in  Figure  3. 


Figure  3.  Transonic  Characteristic  Curves 
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III,  D,  Transonic  Solution  (cont,  ) 


The  derivation  of  the  characteristic  equations  for  the  determination  of  the  particle 
properties  along  the  streamlines  is  shown  in  Appendix  C,  where  is  the  particle - 
stream  function. 


p.  =  _E_ 

ox  u 


u  du  =A(u  -u  )  dx 
P  P  g  P 


V  dv  =A(v  -V  )dy 
P  P  g  P 


u  dh  =-^AC  (T^-T  )dx 
p  P  3  P  g 


d'if  =  0 
p 


This  system  of  equations  can  be  solved  by  assuming  that  all  variable  coefficients 
are  average  values  and  by  putting  the  equations  into  the  following  finite  difference 
form  suitable  for  computer  evaluation: 

Ax  u 
P 


Au  =  A 
P 


Av  =  A 
P 


u  -  u  \ 


Ah  = 
P 


4Acfl 
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III,  D,  Transonic  Solution  (cont.  ) 

may  be  evaluated  along  the  sink  flow  line  to  establish  a  value  of  which, 
from  equation  (68),  remains  constant  along  the  streamline,  thus  determining 
l^p  along  the  supersonic  starting  line.  By  using  the  procedure  shown  in  Appendix 
C,  ■'j^p  along  the  sink  flow  line  may  be  evaluated  as  follows. 

dij^  =  )  dx  +  )  dy  (7  3 

^p  p  X  p  y 


=  d  P  /  u  dy  -  V  dx 
P  P  P  P 


The  value  of  p  along  the  sink  flow  line  may  be  found  by  using  the  sink-flow  solu¬ 
tion,  and  Up  and  Vp  may  be  determined  by 

u  =  Ku  (75) 

P  g 

V  =  Kv  (76) 

P  g 


where  u  and  v  are  determined  from  the  transonic  approximation  along  the  sink 
g  g 

flow  line,  and  K  is  derived  from  the  sink-flow  solution  along  the  sink  flow  line. 
The  value  of  Pp  along  the  supersonic  starting  line  may  be  found  by  solving  (74) 
for  Pp. 


**P  O  (u  dy  -  V  3j^ 

P  p 


The  values  of  P  and  p  along  the  supersonic  starting  line  can  be  found  by 

O  O 

numerically  integrating  the  two  characteristic  equations  which  are  valid  along  the 
gas  streamlines  from  the  sink  flow  line  to  the  starting  line.  This  will  establish 
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III,  D,  Transonic  Solution  (cont. ) 


a  system  of  gas  streamlines  in  the  transonic  region  in  the  same  manner  as  the 
particle  streamlines  illustrated  in  Figure  3  are  established.  By  using  the  deriva¬ 
tion  of  the  characteristic  equations  in  Appendix  C,  the  equations  necessary  to  de¬ 
termine  P  and  p  are 
g 


^  -  _Ii 


p  fu  du  +  V  dv  ]+  dP  =  -  Ap  r  (u  -  u  )  dx  +  {v 

^gL  ggggJ  g  PLgP  g 


Vp)  dyj 


Ap  Bdx 

__2 _ 

P  u 
g  g 


(78) 

(79) 

(80) 


For  the  gas  in  transonic  flow,  the  equations  in  finite-difference  form  are 
.  V 

^  =  -S- 
Ax  — 


=  -A'T  r (a  -  u  )  Ax  +  (v  -  V  )  Ay  -  p  Tu  Au  +  v  Av  I 
g  ^pL'  g  p'  g  P  '^gL  g  g  g  gJ 


Ap  Bp  Ax 
Ap  =  I  AP - E - 8 - 

g  zrr  g 

a 


Y  P  u 
g  g  g 


(81) 

(32) 

(8?) 


Thus  Pg  and  p^  can  be  determined  along  the  starting  line  to  the  same  order  of 

approximation  as  was  involved  in  finding  u  and  v  . 

g  g 


Once  all  the  flow  properties  have  been  found  on  the  starting  line,  the  total 
mass  flow  across  the  starting  line  may  be  found  by  integration.  If  this  total  mass 
flow  is  different  from  the  assumed  subsonic  total  mass  flow,  the  subsonic  total 
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III,  D,  Transonic  Solution  {cont. ) 

mass  flow  should  be  set  equal  to  the  total  mass  flow  through  the  starting  line, 
and  the  entire  subsonic  and  transonic  solution  should  be  repeated.  When  the 
total  mass  flow  through  the  starting  line  finally  equals  the  subsonic  total  mass 
flow,  the  supersonic  starting  line  is  sufficiently  determined  to  continue  with 
the  supersonic  solution  by  the  characteristic  method. 

E.  SUPERSONIC  SOLUTION 


In  supersonic  flow,  the  equations  are  totally  hyperbolic,  and  the 
characteristic  curves  are  all  real.  Hence,  the  flow  field  can  be  determined  by 
a  numerical  technique  based  on  the  characteristic  equations  derived  earlier. 

As  discussed  in  Appendix  B,  to  obtain  the  solution  in  a  region  of  a  flow  field,  an 
initial  data  line  that  is  nowhere  characteristic  and  along  which  all  the  flow  pro¬ 
perties  are  known  must  be  determined.  This  line  is  the  supersonic  starting 
line  discussed  in  section  III,  D  and  must  be  determined  in  a  region  where  the 
Mach  number  is  greater  than  1. 


Once  the  starting  line  has  been  determined,  the  supersonic  solution 
can  be  initiated.  There  are  several  types  of  flow  regions  within  the  supersonic 
regime,  and  each  must  be  handled  differently.  Points  near  the  axis  must  be 
treated  in  a  special  manner  because  of  the  term  (1/y)  in  several  of  the  equations. 
Solid  and  free  boundaries  must  satisfy  the  added  condition  of  known  flow  angle 
and  known  pressure,  respectively.  In  the  region  between  the  limiting  particle 
streamline  and  the  nozzle  wall,  no  particles  are  present  and  the  gas  flow  re¬ 
duces  to  the  case  of  rotational  perfect-gas  flow.  The  flow  near  a  limiting 
particle  streamline  must  also  be  handled  differently,  since  =  0  on  one  side  of 
this  streamline  and  is  finite  on  the  other  side.  The  final  type  of  flow  is  that  of  a 
general  interior  point  where  none  of  the  above  special  situations  prevail.  The 
procedure  for  determining  the  flow  in  the  neighborhood  of  an  interior  point  is 
developed  in  detail  in  this  study.  The  types  of  points  discussed  are  illustrated 
in  Figure  4. 
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Figure  4.  Special  Cases  of  Supersonic  Flow 
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III,  E,  Supersonic  Solution  (cont.  ) 

The  nozzle  wall  is  determined  by  specifying  x  and  y  coordinates  and 
the  slope  at  a  discrete  number  of  points,  as  illustrated  in  Figure  4,  A  left  Mach 
line  is  then  forced  to  go  through  each  wall  point  to  determine  the  flow  properties 
at  the  wall  point,  as  illustrated  at  point  3  for  a  solid  boundary  with  particles  and 
at  point  7  for  a  solid  boundary  with  no  particles.  A  right  Mach  line  is  then  pro¬ 
jected  out  into  the  flow  field  from  each  wall  point  to  continue  the  solution.  The 
size  of  the  characteristic  net  can  be  regulated  by  varying  the  spacing  of  the 
specified  wall  points. 

The  starting  line  is  broken  up  into  several  discrete  points  from  which 
the  characteristic  solution  is  initiated.  The  size  of  the  intervals  chosen  along  the 
starting  line  determines  the  size  of  the  characteristic  net  near  the  nozzle  throat 
where  property  gradients  are  generally  large.  A  relatively  small  interval  is  re¬ 
quired  along  the  starting  line. 

A  characteristic  coordinate  system  is  indicated  in  Figure  4  for  the 
solution  as  it  proceeds  down  the  nozzle.  As  indicated,  right  Mach  lines  are 
numbered  from  1  through  i,  and  left  Mach  lines  are  numbered  from  1  through  j. 
Right  Mach  line  1  initiates  on  the  starting  line  at  the  nozzle  axis,  and  a  right 
Mach  line  of  increasing  number  initiates  from  each  point  on  the  starting  line  up 
to  the  wail,  where  a  right  Mach  line  then  initiates  from  each  wall  point.  Left 
Mach  line  1  initiates  on  the  starting  line  at  the  wall,  and  a  left  Mach  line  of  in¬ 
creasing  number  passes  through  each  point  on  the  starting  line  until  the  nozzle 
axis  is  reached.  Left  Mach  lines  then  are  initiated  at  each  point  on  the  nozzle 
axis  where  a  right  Mach  line  from  a  nozzle  wall  point  intersects  the  axis,  as 
indicated  in  Figure  4.  Thus,  every  point  in  the  characteristic  solution  can  be 
specified  by  a  coordinate  (i,  j)  in  the  characteristic  coordinate  system. 

Now  that  the  characteristic  nets  in  the  supersonic  region  are 
specified,  it  is  necessary  to  develop  a  numerical  solution  procedure  for  each  of 
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III,  E,  Supersonic  Solution  (cont.  ) 

the  points  discussed  above  and  an  overall  calculation  scheme  for  determining 
the  flow  field  in  the  entire  nozzle.  Point  1,  the  general  interior  point,  has  been 
considered  in  great  detail,  and  the  finite-difference  equations  necessary  to 
solve  for  such  a  point  are  discussed  in  section  III,  F. 

F.  FINITE-DIFFERENCE  EQUATIONS 

As  discussed  in  Appendix  C,  it  was  necessary  to  solve  the  particle 
continuity  equation  by  a  mass  balance  on  the  characteristic  control  volume.  The 
remaining  equations  were  found  to  constitute  a  characteristic  system  that  is 
hyperbolic.  The  finite-difference  form  of  these  characteristic  equations  was 
derived  in  Appendix  G  for  the  characteristic  net  shown  in  Figure  5.  Points  3,  4, 
and  5  were  located  by  geometric  considerations  of  the  characteristic  curves,  as 
discussed  in  Appendix  H. 
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III,  F,  Finite -Difference  Equations  (cont,  ) 

To  account  for  more  than  one  size  or  species  of  particle,  six  discrete 
particles  were  accounted  for  by  introducing  the  streamline  of  each  additional 
particle  as  a  characteristic  curve  along  which  the  particle  compatibility  equations 
are  valid.  Particle  parameters  appearing  in  gas  compatibility  equations  were 
summed  for  the  six  particles  to  account  for  their  effects  on  these  equations. 

The  equations  necessary  to  calculate  the  gas  and  particle  properties 
at  point  4  were  derived  as  shown  in  Appendix  G  and  are  discussed  in  detail  for  a 
general  interior  point  in  Appendix  H.  These  equations  include  the  properties  at 
point  4,  which,  as  a  first  approximation,  were  assumed  to  be  the  same  as  the 
properties  at  the  intersection  of  the  initial  data  line  and  the  characteristic  line 
along  which  the  particular  parameter  was  evaluated.  For  subsequent  calculations, 
the  properties  at  point  4  were  considered  to  be  equal  to  the  results  of  the  previous 
calculations  at  point  4.  In  this  manner,  an  iteration  procedure  was  developed  for 
the  flow  properties  at  point  4.  A  numerical  example  also  was  derived,  as  shown 
in  Appendix  I,  to  demonstrate  the  procedure  for  determining  a  general  interior 
point. 


G.  INTERIOR  POINT  NUMERICAL  EXAMPLE 

To  check  the  validity  of  the  theory  and  the  correctness  of  the  IBM 
7090  computer  program,  a  numerical  example  was  calculated  for  a  general  in¬ 
terior  point  (Appendix  I).  A  solid  propellant  containing  aluminum  was  selected 
as  a  sample  for  this  calculation.  The  particles  in  the  exhaust  stream  contained 
aluminum  oxide  particles  of  six  different  sizes.  The  calculation  was  made  for 
a  point  1  ft  downstream  of  the  nozzle  throat  and  0.  5  ft  from  the  nozzle  axis  in 
an  axisymmetric  nozzle. 

Tabulated  values  of  f  and  g  as  functions  of  Reynolds  number  are 
valid  for  spheres  of  any  chemical  species  and  size  and  may  be  considered  to  be 
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III,  G,  Interior  Point  Numerical  Example  (cont.  ) 

the  proper  relationship  for  future  calculations.  The  values  of  Y  c  as 

functions  of  gas  temperature  were  estimated  for  this  particular  gas  and  should 
be  re-evaluated  for  future  cases.  The  temperature -enthalpy  relationship  for 
the  aluminum  oxide  shown  in  Appendix  I  may  be  used  for  future  calculations  in¬ 
volving  aluminum  oxide  particles.  The  remaining  flow  properties  were  estimated 
to  allow  a  numerical  calculation  to  be  made  and  are  not  indicative  of  the  type  of 
flow  field  that  would  normally  occur. 

The  numerical-solution  procedure  shown  in  Appendix  H  was  followed 
in  achieving  both  the  calculated  and  the  IBM  7090  computer  solutions.  The  re¬ 
sults  obtained  through  the  two  methods  were  in  agreement,  indicating  that  the 
computer  solution  is  correctly  programed.  The  flow  properties  calculated  at 
the  new  point  were  all  found  to  be  of  the  correct  order  of  magnitude  and  chang¬ 
ing  in  the  proper  direction  for  flow  down  a  diverging  supersonic  nozzle.  There¬ 
fore,  the  theoretical  development  appears  accurate  and  the  computer  program 
appears  to  correctly  calculate  the  flow  field  at  a  general  interior  point. 

IV.  CONCLUSIONS 


The  model  developed  in  this  Technical  note  for  the  two-dimensional  and 
axisymmetric  flow  of  gas -particle  mixtures  appears  to  adequately  represent 
the  actual  flow  field  for  spherical  particles  inthe  exhaust  gases.  Such  is  the 
case  when  the  particles  are  aluminum  oxide.  However,  not  all  particles  are 
necessarily  spherical,  or  even  approach  it,  and  care  must  be  exercised  when 
considering  such  flows. 

The  theory  developed  here,  when  completed,  will  enable  the  determina¬ 
tion  of  the  performance  characteristics  of  conventional  de  Laval  nozzles  and 
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IV,  Conclusions  (cont.  ) 

unconventional  nozzles,  such  as  the  plug  nozzle.  Particle  trajectories  and  im¬ 
pingement  points  can  be  accurately  determined.  Nozzle  thrust  and  specific  im¬ 
pulse  can  be  calculated,  thus  allowing  nozzle  evaluation  studies  to  be  made  to 
determine  how  to  design  nozzles  for  gas -particle  mixtures. 

V.  RECOMMENDATIONS  FOR  FUTURE  WORK 


The  work  described  herein  is  the  result  of  a  study  of  the  flow  of  gas- 
particle  mixtures  in  two-dimensional  and  asixymmetric  rocket  nozzles.  The 
program  is  far  from  complete,  however^  Further  work  needs  to  be  programed 
for  the  following  specific  items: 

1.  A  more  detailed  analysis  of  the  subsonic  solution, 

2.  A  more  detailed  analysis  of  the  transonic  solution, 

3.  Development  of  a  solution  procedure  for  the  special  cases  of  supersonic 
flow  mentioned  in  section  III,  E, 

4.  Development  of  a  computer  program  for  subsonic,  transonic,  and  super¬ 
sonic  flow,  uniting  into  one  program  the  solutions  relating  to  the  entire 
nozzle, 

5.  Analysis  of  particle  trajectories, 

6.  Development  of  an  analysis  procedure  for  unconventional  nozzles, 

7.  Consideration  of  heat  transfer  to  the  wall, 

8.  Consideration  of  chemical  nonequilibrium  in  the  gas  phase, 

9.  A  systematic  computer  analysis  of  the  effects  of  nozzle  geometry, 

10.  A  systematic  computer  analysis  of  the  effects  of  gas-dynamic  parameters, 

11.  Correlation  with  an  experimental  program. 

12.  Re-evaluation  of  the  theoretical  analysis  based  on  the  experimental  com¬ 
parison. 

13.  Development  of  a  complete  method  for  designing  any  type  of  nozzle  for 
gas -particle  systems. 
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VI.  NOMENCLATURE 

A.  ENGLISH  SYMBOLS 

A  Particle  parameter  in  characteristic  equations 
Sg  Parameter  in  particle  finite-difference  equations 
area  Area  of  one -dimensional  gas -particle  system  nozzle 
a  Speed  of  sound  in  a  perfect  gas 

a*  Speed  of  sound  at  sonic  conditions 

a_  Coefficient  of  system  of  quasi -linear  partial -differential  equations 
B  Parameter  in  gas  characteristic  equations 

Parameter  in  particle  finite -difference  equations 

Coefficient  of  system  of  quasi -linear  partial -differential  equations 
C  Parameter  in  characteristic  equations 


C^  Drag  coefficient 
D 

Parameter  in  particle  finite -difference  equations 
C{<y)  Characteristic  curve  parametric  representation 


Coefficient  of  system  of  quasi -linear  partial -differential  equations 

Specific  heat  at  constant  pressure  of  the  gas 

Spe  cific  heat  at  constant  pressure  of  the  particles 

Constant -lag  one -dimensional  flow  parameter 


Vector  drag  force 

Dq  Parameter  in  gas  finite -difference  equations 

Parameter  in  characteristic  compatibility  equations  derivation 
ds  Differential  arc  length  in  radius  of  curvature  calculation 


dx  Differential  distance  along  nozzle  axis 
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VI,  A,  English  Symbols  (cont.  ) 


dy 
d(  ) 

E 

T 

f 


fo(x) 

f2(x) 

g 

h 

h 

g 

h 

P 

K 


Ki,  Kjj 


k 

g 

k. 


L 


L. 

1 


KT 


Mj.  Mjj 


Differential  distance  normal  to  nozzle  axis 
Total  differential  of  a  quantity 
Constant-lag  one -dimensional  flow  parameter 
Vector  force 

Ratio  of  to  Cj^  for  Stokes'  flow 

Coefficient  in  power  series  solution  for  ^  (x,  y)  in  transonic  flow 

Coefficient  in  power  series  solution  for  ^  {x,  y)  in  transonic  flow 

Coefficient  in  power  series  solution  for  ^  {x,  y)  in  transonic  flow 

Ratio  of  Nu  to  Nu  for  Stokes'  flow 

Particle  heat -transfer  film  coefficient 

Enthalpy  of  the  gas 

Enthalpy  of  the  particles 

One -dimensional  particle  velocity  lag 

Parameters  in  gas  finite -difference  equations 

Thermal  conductivity  of  the  gas 

Parameter  in  characteristic  compatibility  equations  derivation 
One -dimensional  particle  thermal  lag 

Notation  for  a  system  of  quasi -linear  partial  -  differential  equations 
Parameter  in  particle  mass  balance 
Parameter  in  particle  mass  balance 
Mach  number 

One -dimensional  constant-lag  modified  Mach  number 
Parameters  in  gas  finite -difference  equations 
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VI,  A,  English  Symbols  (cont.  ) 


m 

m 

P 

Nu 

n 

n 

n. 

1 

P 

g 

P. 

1 

Pr 


Q 

Q 


I 


.  Q 


II 


q 

R 


Re 


^II 


S  ,  S 
I’  II 


Slope  of  a  line  in  geometrical  characteristic  net 
Density  of  a  particle  per  unit  volume  of  particle 
Nussult  number 

Index  in  system  of  quasi -linear  partial -differential  equations 

Unit  normal  vector 

Unit  vector  in  ith  direction 

Gas  pressure 

Any  property  in  finite -difference  equations 
Prandtl  Number 

Heat  added  to  gas -particle  system  control  volume 
Parameters  in  gas  finite -difference  equations 
Velocity  vector 
Gas  constant 
Reynolds  number 

Parameters  in  gas  finite -difference  equations 
Particle  radius 

Parameters  in  gas  finite -difference  equations 
Throat  wall  location  in  transonic  flow  model 
Parameter  in  gas  finite -difference  equations 
Gas  temperature 

Gas  stagnation  chamber  temperature  in  one -dimensional  model 
Particle  temperature 

Particle  stagnation  chamber  temperature  in  one -dimensional  model 
Time 
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VI,  A,  English  Symbols  (cont.  ) 


U 


u 


u 

p 

u  . 
pi 

uJ 

11  j 


n+1 


qJ 

t 

V 

V 

V 

V 


w 


w 


Y 

y 


Parameter  in  gas  finite -difference  equations,  and  nondimensional 
x-direction  velocity  component  in  transonic  flow 

Perturbation  x-direction  velocity  component  in  transonic  flow 

x-direction  velocity  component  in  transonic  flow 

x-direction  gas  velocity  component 

Gas  velocity  component  in  ith  direction 

x-direction  particle  velocity  component 

Particle  velocity  component  in  ith  direction 

Flow  property  in  system  of  quasi -linear  partial-differential  equations 
Value  of  u^  after  n  +  1th  iteration 
Value  of  ul  after  nth  iteration 

Parameter  in  gas  finite -difference  equations  and  nondimensional 
y -direction  velocity  component  in  transonic  flow 

Volume  of  gas -particle  system  control  volume 

Perturbation  y -direction  velocity  component  in  transonic  flow 

y -direction  velocity  component  in  transonic  flow 

y -direction  gas -velocity  component 

y-direction  particle -velocity  component 

One -dimensional  gas  mass -flow  rate 

One -dimensional  particle  mass -flow  rate 

Parameter  in  gas  finite -difference  equations 

Coordinate  along  nozzle  axis 

Direction  cosine  of  characteristic  curve 

Parameter  in  gas  finite -difference  equations 

Coordinate  normal  to  nozzle  axis 
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A,  English  Symbols  (cont.  ) 


Yg  Radius  of  nozzle  throat  in  transonic  flow  model 

y  Direction  cosine  of  characteristic  curve 

Z  Parameter  in  gas  finite -difference  equations 

z  Coordinate  normal  to  x-y  plane  in  nozzle 


B. 

a 


g 


’g 

A(  ) 


Ax 

Ay 

b 


GREEK  SYMBOLS 

Mach  angle,  gas  viscosity  parameter,  or  Sauer -line  velocity 
gradient 

Initial  data  curve  in  method  of  characteristics 
Solution  curve  in  method  of  characteristics 
Specific  heat  ratio  of  the  gas 

One -dimensional  constant-lag  modified  specific  heat  ratio 
Finite  difference  of  a  quantity 
Finite  change  along  nozzle  axis 
Finite  change  normal  to  nozzle  axis 

Denotes  two-dimensional  or  axisymmetric  flow  in  Sauer's 
approximation 

Location  of  coordinate  system  in  transonic  flow 
Slope  of  characteristic  curve  for  two  variables 
Slope  of  characteristic  curve  for  two  variables 
Slope  of  characteristic  curve  for  n  variables 
Angle  of  the  tangent  to  a  gas  streamline 
Angle  of  the  tangent  to  a  particle  streamline 
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VI,  B,  Greek  Symbols  (cont.  ) 


\ 


TT 


I 


Linear  multiplying  factors  in  method  of  characteristics 
Angle  of  a  line  in  geometrical  characteristics  net 
Gas  viscosity 
Gas  viscosity  parameter 
3.  14159. . . 

Gas  density  per  unit  volume  of  gas 
Particle  density  per  unit  volume  of  gas 
Radius  of  curvature  in  transonic  flow 

Indicates  a  summation 


O’  Denotes  type  of  flow,  two  dimensional  or  axisymmetric 

Oy,  Denotes  type  of  flow,  two  dimensional  or  axisymmetric 

T  Angle  for  determining  radius  of  curvature 

^(x,  y)  Gas  velocity  potential  function  in  transonic  flow 

y  Particle  stream  function 

P 

C.  SUBSCRIPTS 

G  Gas  streamline 

g  Gas  property 

go  Gas  property  at  stagnation  chamber  conditions 

i  Index 

ij  Index 

j  Index 

m  Midpoint  of  initial  data  line 
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VI,  C,  Subscripts  (cont.  ) 

p  Particle  property 

po  Particle  property  at  stagnation  chamber  conditions 

S  Particle  streamline 

Stokes'  Stokes'  flow  regime  where  Rg  ^  ^ 

t  Nozzle  throat 

I  Right  running  Mach  line 

II  Left  running  Mach  line 

D.  SUPERSCRIPTS 

j  Index  for  flow  variables 

♦  Sonic  conditions 

Average  value  of  a  quantity,  or  a  vector  quantity 

Velocity  components  in  transonic  flow  model 

E.  OTHER 
Vector  operator 

Total  derivative  with  respect  to  time 
Partial  derivative  with  respect  to  time 
Slope  of  a  line  in  x-y  plane 
Natural  logarithm 
Partial  derivative  with  respect  to  x 
Partial  derivative  with  respect  to  y 

Relative  velocity  vector  between  gas  and  particle~(qg  - 
Absolute  value  of  a  quantity 


[CTj 
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VI,  Nomenclature  (cont.  ) 


g 

Pr 

r 


P 

R 


y 

a 


UNITS  OF  PHYSICAL  PARAMETERS 

Specific  heat  at  constant  pressure 
Ratio  of  to  for  Stokes'  flow 
Ratio  of  Nu  to  Nu  for  Stokes'  flow 
Particle  enthalpy 

Particle  density  per  unit  volume  of  particle 
Gas  static  pressure 
Prandtl  number 

Particle  radius,  based  on  spherical  particles 

Gas  constant 

Reynolds  number 

Gas  static  temperature 

Particle  static  temperature 

Axial  gas 'Velocity  component 

Axial  particle -velocity  component 

Normal  or  radial  gas -velocity  component 

Normal  or  radial  particle -velocity  component 

Axial  coordinate 

Normal  coordinate 

Exponent  of  gas -viscosity  expression 
Specific  heat  ratio  of  the  gas 
Gas  viscosity 

Coefficient  of  gas -viscosity  expression 


Btu/lbm-®R 
dimensionle  ss 
dimensionless 
Btu/lbm 
Ibm/cu  ft 
Ibf/sq  in. 
dimensionless 
microns 
Btu/lbm. ‘R 
dimensionless 
•R 
•R 

ft/sec 
ft/ sec 
ft/sec 
ft/ sec 
ft 
ft 

dimensionle  s  s 
dimensionless 
Ibm/ft-sec 
Ibm/ft-se': -(®R)° 
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VI,  F,  Units  of  Physical  Parameters  (cont. ) 


Gas  density 

p  Particle  density  per  unit  volume  of  gas 
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DERIVATION  OF  EQUATIONS- GOVERNING  GAS-PARTICLE  MIXTURES 


A  gas -particle  mixture  is  described  by  the  equations  for  conservation  of  mass, 
conservation  of  energy,  momentum  balance,  equation  of  state  for  the  gas, 
particle  drag,  particle  heat  balance,  and  particle  enthalpy-temperature  relation¬ 
ships.  This  appendix  derives  these  equations  on  the  basis  of  following  assump¬ 
tions: 

1.  Only  one  particle  size  is  present,  or  some  type  of  integrated  particle 
distribution  is  performed,  and  the  particles  are  spherical. 

2.  The  total  mass  of  the  gas -particle  mixture  is  constant. 

3.  The  total  energy  of  the  gas -particle  mixture  is  constant. 

4.  The  internal  temperature  of  the  particles  is  uniform,  and  the 
particle  specific  heat  is  constant, 

5.  The  gas  and  particles  exchange  thermal  energy  by  convection  only, 

6.  The  gas  obeys  the  perfect-gas  law,  has  a  constant  molecular  weight, 

and  constant  specific  heats. 

7.  All  external  forces  except  pressure  of  the  gas  and  drag  of  the  particles 
are  neglected. 

8.  The  gas  is  inviscid  except  for  the  drag  it  exerts  on  the  particles. 

9.  The  particles  do  not  interact  with  each  other. 

10,  The  volume  occupied  by  the  particles  is  negligible. 


The  following  derivations  of  the  equations  are  based  on  studies  by  Kliegel:’!' 


♦  Kliegel,  J.  R.  and  G.  R,  Nickerson,  Flow  of  Gas -Particle  Mixtures  in  Axially 
Symmetric  Nozz^s,  STL  TM-7106-00i3-MU-000;  also  ARS  preprint  1713-bl, 
presented  at  the  Propellants  Combustion  and  Rockets  Conference,  April  26-28, 
1961. 
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To  begin  the  derivation,  a  control  volume  V  enclosed  by  a  surface  A  is 
considered.  This  control  volume  is  illustrated  by  Figure  A-1  below. 


The  control  volume  is  fixed  in  space  and  the  gas-particle  mixture  flows  through 
it.  The  mass  flux  through  the  surface  A  is 


q  •n  + 
g 


Pp  ‘Ip 


^n 


dA 


(A-i) 


where  is  the  mass  of  particles  per  unit  volume  of  gas. 
change  of  mass  in  V  is 


The  time  rate  of 


A  mass  balance  for  V  results  in 


(A-E) 


(A -3) 
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When  the  divergence  theorem  is  used,  the  surface  integral  of  Equation 
A-3  can  be  transformed  to  a  volume  integral,  and  Equation  A-3  becomes 


r 

J  .8‘  at 


+  div  P  q  +  P  q 
g  g  P 


dV  =  0 


(A -4) 


Since  this  equation  must  be  true  for  any  size  volume,  it  can  be  written  as 


(p^qjl  t[^+div(ppqp)  = 


However,  there  is  no  mass  interchange  between  the  gas  and  particles  in  V; 
consequently,  each  member  of  Equation  A- 5  must  be  identically  zero.  Thus, 
the  conservation  of  mass  equation  for  V  take  the  form 


arr  ^  Pa  ° 


{A*  6i 


j^+div  p  q  =  0 


The  momentum  flux  through  A  in  the  ith  direction  is 


f  [K  v)  ’g*"^  (  Vpi)  ’p’". 


{A -7^ 


(A -8) 


The  time  rate  of  change  of  momentum  in  V  in  the  ith  direction  is 


/[-^(vgi)  (vpi)] 


(A -9) 
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Through  the  assumption  that  all  body  forces  except  the  pressure  of  the  gas 
and  the  drag  are  negligible,  there  are  no  viscous  or  body  forces  acting  on 
surface  A,  and  the  only  forces  present  are  pressure  forces.  The  component 
of  this  pressure  force  in  the  ith  direction  is 


fp  n.c 

J  g  1 


(A- 10) 


A  momentum  balance  for  V  gives 


V"  ^  _ 

2_,  IP  =  HF  Momentum 


^  -  -1 

"pi  I  S* " 


V""  "p 


+'5t“ip  u ,  n  dv 


'  p  u  . 
p  pi 


(A- 11) 


When  the  divergence  theorem  is  applied  to  Equation  A- 11,  the  expression  for 
the  ith  component  of  the  momentum  balance  becomes 


/  I  Qt  I  p  u  .  +  p  u  .  I  +  div  I  p  u  .  q  +  p  u  .  q_ 

J  I  \  g  gi  P  pi  I  ^  g  gi  g  P  pi  T 


(A- 12) 


Op  1 

Dif-J 


dV  =  0 


Since  this  equation  must  hold  for  any  volume,  the  momentum  balance  becomes 


fp  u  .  +  p  u  .  ]  +  div  fpu.q  +  pu.q\  +  ?  =  0 

ft  r  g  gi  p  pi  I  I  ^g  gi  g  p  pi  P I  ^ 


(A- 13) 
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The  flux  of  enthalpy  and  kinetic  energy  through  A  is 


the  time  rate  of  change  of  enthalpy  and  kinetic  energy  in  V  is 


h  + 

g 


1 


r 


d 

+  w~ 


dV 


The  time  rate  of  heat  addition  to  V  is 


An  energy  balance  for  V  gives 


When  the  divergence  theorem  is  applied  to  Equation  A -17 


(A- 14'? 


(A-15) 


(A- 16) 


{A- 17) 


(A-IS) 
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Since  this  equation  must  hold  for  any  volume,  the  energy  balance  becomes 


(A-19) 


With  the  assumption  that  the  gas  obeys  the  perfect  gas  law  and  has  a 
constant  molecular  weight  and  constant  specific  heat,  the  equation  of  state 
of  the  gas  can  be  written  as 


P  =  p  R  T  (A-20) 

g  g  g 

Next,  a  spherical  particle  of  constant  size  is  considered  in  the  gas 
flow  field  as  shown  in  Figure  A-2. 


Figure  A-2.  Particle  Momentum  and  Heat  Transfer  Model 

Remembering  the  assumption  that  all  body  forces  except  for  pressure  of 
the  gas  and  particle  drag  are  negligible,  the  only  force  exerted  on  the  particle 
is  a  drag  force,  13|,  due  to  the  relative  motion  between  the  gas  and  particle. 

A  particle  momentum  balance  gives 

°£'-af  ("'pT ’■■ps) 
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An  energy  balance  for  the  particle  gives 


at 


=  -  h4Trr^  f  T  -  T  ) 
P  '  P  g' 


k  kg 

h  =  ^  Nu  =  -SI 
2r  r 

P  P 


Q  =  m  iirr^  h 
P  P  3  p  p 


dh  n  g 

-JE.=  -  3 

dt 


m  r 
P  P 


T  k 


T  ) 
g  / 


{A-29) 

(A-30) 


{A-31) 

(A-32) 


The  temperature  of  the  particle  depends  on  its  state;  either  a  liquid,  a  liquid 
in  the  process  of  solidifying,  or  a  solid,  and  is  uniquely  related  to  the  enthalpy, 

Tp  =  £  (  hp)  (A.33) 


The  gas -particle  system  is  now  completely  defined  by  Equations  A-6, 

A-7,  A-13,  A-19,  A-20,  A-28,  A-32,  and  A-33.  For  adiabatic  steady  flow,  the 
two-dimensional  and  axisymmetric  forms  of  these  equations  are  obtained  by  drop¬ 
ping  all  time  derivatives  and  expanding  the  above  system  of  equations.  By  let- 
tinge;  =  1  for  two-dimensional  flow,  O  -  y  for  axisymmetric  flow,  x  denote  the 
longitudinal  axis,  and  y  denote  the  vertical  or  radial  axis,  these  equations  re¬ 
duce  to  the  following; 


p(u)  +p(v)+u(p)+v(p)  =-a  ^SL^a  (A-34) 

g  X  g  g  y  g  ^g  X  g  ^g  y  y  ' 

p  (u)  +  p  (v)  +  u(p)  +  v  (p)  =-<y  (A-35) 

P  X  '  p'y  p'^p'x  p  '»'pV  y 
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Expanding  Equation  A- 13  and  combining  with  the  expanded  form  of  . 
Equation  A-28  gives 

p  pu  (u  )  +  V  (vi  )  ]+-Ap  (u  -  u  )  +  (P  )  =0 

L.  g'  g'x  g  g'yj/  '  g  P  g'x 

p  Pu{v)  +v  (v)l  +  Ap  (v  -  V  )  +  (P  )  =0 

^g  L  g'  g'x  g  g'yj  ^p  '  g  p'  g'y 


P  P 


The  particle  momentum  and  energy  equations  become 


(A.36) 

{A-37) 

(A~38) 


(u  ) 

+  V 

(u  ) 

*  A  (u  -  u_) 

P  X 

p 

p  y 

g  p 

(v  ) 

+  V 

a  A  (v  -  V  ) 

p  X 

p 

p  y 

g  p 

+  V 

a  -  2  AC  (T  -  T  ) 

p  X 

p 

p  y 

y  P  8 

{A~39) 

!;a-40) 

(A.41) 

('A-42) 


Expanding  the  system  energy  equation  (Equation  19)  and  combining  with 
the  gas  and  particle  continuity  equations,  the  gas  equation  of  state,  the 
expanded  form  of  the  system  monfentum  equation,  and  the  following  gas 
enthalpy -temperature  relationship  gives 


h  =  c  T  +  Constant 

g  Pg  g 


u  (P  )  +  V  (P  )  -  a^  P  u  (p  )  +  V  (p  )  "1  -  Ap  B  =  0 

g  g  X  g  '  g'y  L  g  '^g'x  g  '•'gV  J 


a  =  y  R  T 
g  g 


(A.  43} 

(A.  44} 
CA-45) 
(A-46) 
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The  gas  equation  of  state  and  particle  enthalpjr-temperature  relationship  are 
the  same. 

The  equations  governing  a  gas -particle  flow  field  are  thus  found  to  be  a 
system  of  eight,  quasi -linear,  nonhomogeneous  partial  differential  equations 
of  the  first  order,  shown  as  Equations  A-34,  A-35,  A-36,  A-37,  A-39,  A-40, 
A-41,  and  A-44.  The  constants  A,  B,  and  C  in  these  equations  are  defined  by 
Equations  A -38,  A -42,  and  A -46.  The  gas  equation  of  state  is  given  by 
Equation  A-20,  and  the  particle  equation  of  state  by  Equation  A-33. 
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METHOD  OF  CHARACTERISTICS 

The  governing  differential  equations  often  reduce  to  quasi-linear  partial- 
differential  equations  of  the  first  order  for  functions  of  two  independent  variables. 
A  quasi-linear  partial  differential  equation  of  the  first  order  is  defined  as  one  that 
is  nonlinear  in  the  dependent  variables  but  linear  in  the  first  partial  derivatives  of 
the  dependent  variables.  Such  a  system  of  n  equations  can  be  written  as;* 


0u‘’ 


eu-* 


L.  =  a..  ^  +  b..  ^  +  c. 
1  .  ij  9x  ij  9y  1 


i  =  1,  .  . ,  n 


(B-1) 


where  the  superscript  j  identifies  a  particular  dependent  variable,  and  the  co¬ 
efficients  a  ...  b_.  and  c^  depend  on  X,  y,  u^ . u^.  When  expanded,  this 


system  of  equations  becomes: 
1  .. 


+  “y  +h2"y  +  •  •  •  * ‘>1„  ' 

h  “  »2l  \  *•••**20  "**'’21  "y  *‘’22“y  ^••■+'>2n“y  *  '2'° 

••  •  ?*•  •• 

•  •  •  0© 

'  »nl  *  *„2  "x  *  •  •  ■  V  *  "nl  "y  +  '>„2  “y  +  ’  •  •  +  '’nn  "y  +  « 


When  such  systems  are  hyperbolic,  the  method  of  characteristics  can  be  used 
to  obtain  the  desired  solution.  A  treatise  on  the  method  of  characteristics  ap¬ 
plied  to  fluid-flow  problems  has  been  written  by  Courant  and  Friedrichs,  ** 
and  this  discussion  is  based  on  that  treatise.  The  simplify  presentation,  the 
theory  will  first  be  developed  for  a  system  of  two  equations,  and  then  the  re¬ 
sults  will  be  applied  to  systems  of  n  equations. 

^  In  accordance  with  accepted  convention,  when  an  index  is  repeated,  the  sum¬ 
mation  is  carried  out  with  respect  to  that  index. 

**  Courant,  R,,  and  K.  O.  Friedrichs,  Supersonic  Flow  and  Shock  Waves,  New 
York,  Interscience  Publishers,  1948. 
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The  concept  of  characteristics  arises  from  the  following  relationships. 

A  linear  combination  af  +  bf  of  the  two  partial  derivatives  of  a  function  f  can 

X.  y 

be  considered  as  the  derivative  of  f  in  a  direction  given  by  dx:dy  =  a;b.  Figure 
B-1  illustrates  this  point.  In  the  figure,  the  curve  C(o’)  is  given  in 


Figure  B-1  Characteristic  Curve 

parametric  representation  as  a  function  of  some  parameter  Cf.  Thus,  at  any 

point  (x,  y),  x^  and  y^  are  the  direction  cosines  of  the  curve  C  {(T).  If  x^^: 

y-;  ia:b,  then  df/do  =  f  (dx/dc7)  +  f  (dy/do)  =  f  X  +f  Y  ,  and  af  +  bf 

7<y  X  y  ^'xCTyff  xy 

is  the  total  derivative  of  f  along  this  curve. 

Next,  a  system  of  two  equations  for  the  two  dependent  variables  u{x,  y)  and 
/(x,  y)  is  considered: 


L  =  A  u  +  B  u  +  C.  v  +  D  y  +  E  =  0 

1  lx  ly  lx  ly  1 

L,  =  A,  u  +  B-  u  +  C,  v  +  D-  v  +  E-  =  0 

2  2x  2x  Zy  2 


(B-3) 
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In  addition  to  these  equations,  there  are  two  additional  equations, 


du  =  u  dx  +  u  dv 
X  y  •' 

(B-4) 

dv  =  V  dx  +  V  dy 
X  Y  ’ 


that  are  valid  for  continuous  functions  having  continuous  derivatives.  The 
equations  of  B-4  are  applicable  to  the  functions  considered  here  because  these 
functions  are  assumed  to  be  continuous  and  have  continuous  derivatives.  A 
linear  combination,  L,  of  and  is  sought: 

L  =  +^2  ^2  '  ° 

so  that,in  this  differential  expression  for  L,  the  partial  derivatives  of  u  and  v 
combine  to  give  total  derivatives  of  u  and  v  in  the  same  direction.  Such  a 
direction,  if  it  exists,  depends  on  the  point  (x,y),  as  well  on  the  values  of  u  and 
V  at  that  point,  and  is  called  characteristic.  Expanding  L  as  indicated  in 
Equation  B-5  results  in  the  following: 

(XiAf  +\2'^2^“x  '*'^2^2^^x 

+  (\jD^  +\2D2)  Vy +  {\jEj  +  \2E2)  =  0  (B-b) 

Suppose  the  characteristic  direction  is  given  by  curve  C  {a)  with  direction 
cosines  x^  and  y^  .  Then,  for  u  and  v  in  the  expression  for  L  to  be  differ¬ 
entiated  in  the  direction  of  C  (a)  only,  by  the  analysis  above  concerning  af^ 

+  bfy,  the  following  proportionality  must  be  satisfied: 

(XlAf  +  \2A2)  :  +\2B2)  =  \^2^  ’  ^^2^2^  “  *tr ’Xff 

since  the  coefficients  of  the  derivatives  u  ,  u  ,  v  •  and  v  in  L  are  given  by  the 

x  y  x  y 

respective  members  of  the  above  proportion. 
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If  the  functions  u  and  v  satisfy  the  differential  equations  and  at 
the  point  (x,y),  the  following  two  homogeneous  linear  equations  for  and 
are  obtained: 


(B-8) 


For  B-8  to  have  a  nontrivial  solution,  the  determinant  of  the  coefficients  of 
X^  and  X2  in  the  above  equations  must  vanish. 


<^1%  <Va  '®2V 


(B-9) 


<^1^2  -  ^z^i^ya  ■  •  ^2^1)]  Vd 


+  (B^D2  -  B^D^)  0 


ay^  “  2bx^y^  +cx^  »  0. 


(B-10) 


(B-ll) 


If  ac-b^>  0,  then  no  real  solutions  for<^  *  ^  exist,  and  the  char- 
acteristic  curves  C  (d)  are  imaginary.  Differential  equations  that  result 
in  imaginary  characteristic  curves  are  called  elliptic.  If  ac-b^  »  0,  one 
real  characteristic  direction  exists  through  each  point,  and  the  system  is 
called  parabolic.  If  ac-b^<  ,  two  real  characteristic  directions  exist  through 
each  point,  and  the  system  is  called  hyperbolic.  From  here  on,  only  hyper¬ 
bolic  systems  are  considered.  Introducing  the  slope 


■;■  — =  P- 

Xff  dx 


(B.12) 
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Equation  B-11  can  be  written  as 
-  ZbC  +  c  =  0 


(B-13) 


which,  for  hyperbolic  systems,  has  two  distinct  real  solutions  t:  +  and  K.-. 
Hence,  the  two  characteristic  curves  C  (a)  satisfy  the  two  ordinary  differential 

equations 


Since  the  roots  +  and  -  are  functions  of  x,  y,  u,  and  v,  the  hyperbolic 
character  of  the  system  depends  on  the  particular  functions  u  (x,y)  and  v 
(x,  y)  under  consideration.  When  a  solution  u  {x,y)  and  v  (x,y)  is  inserted 
into  Equation  B-14,  the  equations  dy/dx  =  C+  (u,  v,  x,  y)  and  dy/dx  =  C- 
(u,  V,  X,  y)  are  two  ordinary  differential  equations  of  the  first  order  that 
define  two  families  of  characteristic  curves,  or  simply  characteristics,  C  + 
and  C-  in  the  (x,y)  plane. 


In  the  case  of  two  dependent  variables  u  and  v,  four  equations  were 

found  relating  u  ,  u  ,  v  ,  and  v  . 

*  X  y  X  y 


=  A.u 

1  X 

+  B.u 

1  y 

+ 

C .  V 

1  X 

+  D  V 

1  y 

+ 

Ej  =  0 

=  A,u 

2  X 

+  B,u 

2  y 

+ 

C,v 

2  X 

+  D,v 

2  y 

+ 

E^  =  0 

{B-15) 


u  dx  +  u  dy  =  du 
X  y 

V  dx  +  V  dy  =  dv 

X  y  ' 

By  using  Cramer's  Rule  for  a  system  of  n  nonhomogeneous  linear  equations  in 
n  unknowns,  the  partial  derivatives  of  u  and  v  may  be  determined.  Cramer's 
Rule  states  that  if  the  determinant  of  the  coefficients  of  the  unknowns,  D,  is 
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D. 

different  from  zero,  each  unknown  has  a  unique  solution  given  by  > 

where  is  the  determinant  obtained  from  D  by  replacing  the  ith  column 
of  h  by  the  nonhomogeneous  terms,  and  k.  is  any  one  of  the  partial  derivatives 
of  u  and  v.  In  the  present  case  of  two  dependent  variables  discussed  above,  D 
is  found  to  be 


^ 

(A^dy  -  Bjdx)(C^dy  -  D^dx)  B^ 

D  a 

■^2  ®2  ^2  ^2 

= 

(A2dy  -  B2dx){C2dy  -  D2dx)  B2 

dx  dy  0  0 

(dydx  -  dydx)  0  dy  0 

0  0  dx  dy 

0  (idxdy  -  dydx)  0  dy 

D 


(Ajdy  -B^dx)  (C^dy  -  D^dx) 
(A2dy  -  B^dx)  (C2dy  -  D2dx) 


(dy)^ 


(B-16) 


By  comparing  this  Solution  for  O  with  Equation  B-9,  it  is  found  that  D  s  0  along 
characteristic  curves.  This  result  is  valid  along  only  the  characteristic  curves. 
Thus,  the  solution  for  k^  =  D^^/D  will  not  exist  along  the  characteristic  curves 
unless  D.  =  0,  thus  placing  k.into  the  indeterminant  form: 


(B-17) 


Solving  the  determinants  =  0  results  in  a  system  of  compatibility  equations 
equal  in  number  to  the  number  of  partial  derivatives  in  the  original  system  of 
equations.  The  number  of  compatibility  equations  that  are  independent  equals 
at  most  the  number  of  dependent  variables  in  the  original  system.  By  observing 
the  nonhomogeneous  term  in  the  equations  of  B-15,  it  can  be  seen  that  these  com¬ 
patibility  equations  will  contain  only  total  derivatives  of  u  and  v.  Thus,  the  system 
of  partial  differential  equations  has  been  replaced  by  a  system  of  total  differential 
equations  valid  along  the  characteristic  curves. 
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So  far,  a  fixed  solution  u  {x,y)  and  v  (x,y)  has  been  assumed.  How¬ 
ever,  the  equations  of  B-15  no  longer  depend  explicitly  on  this  solution,  since 
all  the  coefficients  are  known  functions  of  x,  y,  u,  and  v.  By  a  slight  change 
in  interpretation,  the  equations  of  B-15  and  the  compatibility  equations  can 
be  considered  as  a  system  of  four  equations  for  the  determination  of  x,  y,  u, 
and  V.  Replacing  the  original  system  (equations  of  B-3)  by  this  character¬ 
istic  system  is  the  basis  of  the  method  of  characteristics. 


The  method  of  characteristics  can  now  be  extended  to  a  system  of  n 
differential  equations  by  analogy  to  the  above  case  of  two  equations.  Equation 
B-1  for  such  a  system  is  repeated  here  as 


As  in  the  case  of  two  equations,  the  following  equations  result  from  the  con¬ 
tinuous  nature  of  the  assumed  solution  for  u'^: 

du-^  «  u^-^dx  +  Uy-^dy  (j  =  1,  . . . ,  n)  (B-18) 

Again,  a  specific  solution  u^  is  assumed,  and  curves  C  ((j)  are  sought  so 
that  a  linear  combination  of  the  differential  equations  can  be  formed 
in  which  differentiations  occur  only  along  the  curves  C  (o). 


L  =  \  L  =  \,L  +  \  L  +  ...  +  \  L  =  0 

11  11  l  z  nn 


(B-19) 


Expanding  L  as  indicated  in  Equation  B-19  yields  the  following: 

(\aii  +\2a2j  +  ..  +  (X^b^^  d-X^b^^  +  ..  +X^.b^j)  u^^ 

+  (Xja^2  +  X2a22  +  •  •  +X^a^2^  ^  ' '  +  ^n^Z^  \  (B-20) 

+  .  . .  +  (Xja^^  +  "x  +  (X^bj^  +  +  . .  +X^b^j^n  ^  ^ 

y 
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Ab  in  the  case  of  two  equations,  the  copfficients  of  and  must  be  in  direct 
proportion  to  the  direction  cosines  and  ^  of  C  (ff)  if  L  is  to  be  differentiated 
only  along  C  (tf).  This  proportionality  takes  the  form 


■*■^2*^21  *  ■ 

+X  b  . )  =  X 
n  nl  0 

=  ^0 

^^2*^22  ^  ■  ‘ 

+\  b  ,)  =  x^ 
n  n2  <J 

•  • 

■  (B-21) 

+\_b,  +  .  . 

2  2n 

•  •• 

•  • 

+X  b  )  =  X  : 
n  nn  0 

Xr 

Equation  B-21,  when  rearranged  withV  as  the  unknown,  takes  the  form 

■  “21^’  ■  '’nl*®*  “  ® 

^  -  '=12*ff>  +^2<^22y®  -  '’22*»>  +’  ’  '  '“n2*ff>  =  “ 


(B-22) 


+''2<“2nyff  '  '=2n^ff>  +•  '  '  ‘’nn*® >  ‘  “ 

which  can  be  written  as 

(a.  .y  -  b.  .*  )  =  0  (j  =  1,,  . .  . ,  n) 

1  iy(J  ij  ff'  ’ 

For  the  solution  of  the  system  of  equations  defined  by  Equation  B-23 
to  be  other  than  zero,  the  determinant  of  the  coefficients  ofV  must 
vanish: 

I  (a..y  -  b...x  )  I  =  0 
'  ij  a  ij  <r  ' 

When  expanded,  the  determinant  takes  the  form 

<»11%  -  Si*®  >  <»12ni  -  '’12*®  >  •••  ‘‘in”®  -  '’ln*®> 

<»2l''®  -  '’21*»>  <>22”®  ■  ’’22*®  >  •  ■  •  '*2n''®  '  '’2n’‘®>  .  o 

t  •  s  • 

•  •  s  • 

•  •  •  • 

•  •  •  Kn^a  "  ‘’nn' a  ) 


(B.23) 


(B-24) 


(B.25) 
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The  expanded  determinant  results  in  an  algebraic  equation  of  the  nth  degree 
for  C  *  IJy  giving  n  roots  which  determine  n  characteristic 

directions.  If  all  n  roots  are  distinct  and  real,  the  system  is  totally  hyper¬ 
bolic.  There  are  then  n  families  of  characteristics  (<y)  satisfying  the 
ordinary  differential  equations 

each  family  covering  the  domain  of  the  (x,  y)  plane  under  consideration. 

As  in  the  case  of  two  equations,  the  following  equations  relating  the  x 
and  y  partial  derivatives  of  u"^  were  found: 


L.  a  a.  .u-^  +  b.  .u-^  +  c.  =  0 

1  ij  x  ij  y  1 

/  i=l,  . . ,  n 

\ j=  i»  • . ,  n 

(From  Equation  B-1) 

u^  dx  +  u^  dy  =  du^ 

X  y 

(j=  1,  . . ,  n) 

(Prom  Equation  B-1  ^ 

Again,  by  using  Cramer's  Rule  for  the  system  of  2n  linear  nonhomogeneous 

equations  in  2n  unknowns,  the  partial  derivatives  of  u'^  can  be  determined. 

By  denoting  any  of  the  partial  derivatives  of  u'^  by  k^,  the  solution  becomes 

k^  •  D^/O,  where  O  is  the  determinant  of  the  coefficients  of  the  equations  of 

B-1  and  B-18.  and  is  the  determinant  obtained  from  D  by  replacing  the 

ith  column  of  D  by  the  nonhomogenous  terms.  Exactly  as  in  the  case  of  two 

equations,  D  is  found  to  be  the  same  as  the  determinant  of  Equation  B-24, 

which  is  equal  to  zero,  and  defines  the  characteristic  curves  C  (ff).  Thus, 

n 

for  kj.  to  exv-'.  must  be  zero,  and  a  system  of  2n  compatibility  equations 
in  terms  of  the  total  derivatives  u'^  are  obtained,  of  which  at  most  n  are 
independent.  Thus  the  general  system  of  n  partial  differential  equations  for 
n  variables  has  been  replaced  by  a  system  of  n  total  differential  equations 
valid  along  the  n  characteristic  curves,  (CT). 

Next,  the  initial  value  problem  is  formulated  for  the  above  system  of 
hyperbolic  differential  equations.  Assume  a  curve  is  given  in  the  (x,y) 
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plane,  e-nd  continuous  values  of  u'^  arbitrarily  prescribed  along Then,  the 
problem  is  to  determine,  in  the  neighborhood  of  P^,  a  solution  uJ  of  the  system 
that  takes  on  the  prescribed  initial  values  along  P^.  By  replacing  the  original 
system  by  the  characteristic  system,  the  system  reduces  to  total  differential 
equations  along  the  characteristic  curves.  In  general,  these  equations  are  non¬ 
linear,  and  each  is  coupled  with  some  of  the  others.  For  this  reason,  a  solu¬ 
tion  by  a  numerical  iteration  technique  becomes  necessary.  The  compatibility 
equations,  each  valid  along  one  or  more  of  the  characteristic  curves,  can  be 
put  in  finite  difference  form,  and  the  equations  of  the  characteristic  curves 
themselves  can  also  be  put  in  this  form.  Then,  by  moving  along  a  characteristic 
curve,  the  initial  values  of  u'^  along  P^  can  be  extended  into  the  domain  enclosed 
by  the  outermost  characteristic  curves  passing  through  the  initial  data  curve 
P^.  By  continuing  in  small  steps  along  the  length  ofP^,  a  new  curve, P^,  can  be 
obtained  with  all  the  values  of  u"^  determined  along  this  curve,  as  shown  in 
Figure  B-2, 


Figure  B-2.  Characteristic  Initial  Value  Line  and  Solution  Net 
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The  above  considerations  result  in  the  concepts  of  domain  of  dependence 
and  range  of  influence.  The  domain  of  dependence  of  a  point  P  is  the  region  in 
the  (x,y)  plane,  bounded  by  the  outermost  characteristics  passing  through  the 
initial  data  line,  in  which  the  solution  of  the  initial  value  problem  can  be 
established.  The  range  of  influence  of  a  point  Q  on  the  initial  value  line  is  the 
totality  of  points  in  the  (x,y)  plane  which  are  influenced  by  the  initial  data  at 
point  Q.  This  region  consists  of  all  points  whose  domain  of  dependence  contain 
the  point  Q;  therefore,  it  is  the  region  between  the  two  outermost  characteristics 
passing  through  the  initial  point  Q.  The  range  of  influence  and  domain  of  de- 


Figure  B-3 

Domain  of  Dependence 


Figure  B-4 
Range  of  Influence 


For  a  solution  to  be  possible,  the  initial  data  line  cannot  be  characteristic 
at  any  place  unless  initial  data  are  given  along  two  intersecting  characteristic 
curves.  Several  types  of  domain  with  different  types  of  initial  data  lines  can  be 
solved:  however,  only  the  domain  in  wnicn  the  initial  data  line  is  nowhere  char¬ 
acteristic  is  of  interest  in  the  present  discussion. 

By  means  of  the  method  of  characteristics  just  described,  many  com¬ 
plicated  systems  of  partial  differential  equations  can  be  solved,  if  the  system 
is  found  to  be  quasi -linear  and  hyperbolic.  These  conditions  are  frequently 
encountered  in  fluid  flow  problems. 
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DETERMINATION  OF  THE  CHARACTERISTICS  OF 
A  GAS -PARTICLE  MIXTURE 

The  complexity  of  the  equations  describing  a  gas -particle  mixture  requires 
that  a  solution  be  obtained  by  other  than  conventional  methods.  The  system  of 
eight  first-order,  quasi-linear,  nonhomogeneous,  partial -differential  equations 
can  generally  be  solved  by  using  the  general  theory  of  the  method  of  characteristics 
as  outlined  in  Appendix  B.  After  an  initial  effort  to  solve  the  entire  system  of 
eight  equations,  it  was  found  that  four  distinct  characteristics  curves  did  exist, 
but  only  seven  distinct  compatibility  relationships  were  found.  Therefore,  the 
particle  continuity  equation  was  solved  separately  by  the  introduction  of  a  particle 
stream  function.  The  resulting  relationship  was  then  considered  to  define  the 
P^’^ticle  density  throughout  the  flow  field,  since  there  was  no  compatibility 
equation  describing  the  variation  of  p^.  The  solution  was  then  found  for  the  particle 
continuity  equation  as  follows; 

(  Vp)  X  *  ?  «'i’pV  y  '  ^  «=-‘> 

The  definition  of  the  stream  function  was  chosen  as 


(^  )  =  ap  u 

p  y  p 

(C-2) 

(i  )  =-ap  V 

P  X  P  P 

(C.3) 

where  <?  =  1  for  two-dimensional  flow,  <T=y  for  axisymmetric  flow,  x  is  the 
longitudinal  axis  direction,  and  y  is  the  vertical  or  radial  direction  for 
two-dimensional  or  axisymmetic  flow,  respectively.  The  definition  of  a  particle 
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streamline  is 


dy  ^  ^  (C-4) 

3x  u 

P 

The  total  derivative  of^  is 

P 

di  =  t  dx  +  dy 
p  X  y  ^ 

di^  =  O'  p  u  dy  -0  p  V  dx 
p  '^p  p  ■'  p  P 

Along  particle  streamlines,  therefore, 

d^  =  0  (C-7) 

P 

The  remaining  seven  equations  form  a  system  of  equations  that  can 
be  solved  by  the  method  of  characteristics.  The  characteristic  curves  are 
found  to  be  the  gas  streamline,  the  particle  streamline,  and  the  two  conven¬ 
tional  gas  Mach  lines.  For  supersonic  flow,  the  system  is  totally  hyperbolic, 
and  for  subsonic  flow  partially  hyperbolic. 

Since  the  seven  quantities  under  consideration  are  assumed  to  be 
continuous  functions,  there  also  exist  seven  equations  of  the  form 

du^  =  u  ^dx  +  u  ^dy  (j  =  1,  .  . ,  7)  (C-8) 

X  y 

where  u^  represents  u,v,p,P  ,u,v,  and  h  .  Thus,  there  exist 

gg  ggPP  P 

the  following  fourteen  equations  for  the  determination  of  the  fourteen 
partial  derivatives  of  the  seven  dependent  variables; 


(C-5) 

(C-6) 
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D  (u  )  +  p  (v  )  +  u  {p  )  +  V  (p  )  =  -Of 

g  X  ^g'  gV  g'^g'x  g'^'g  y  y 

(C-9) 

^g  g'  g'x  g'  g'y  g'x  '^p  '  g  p' 

(C-IO) 

pu(v)  +  pv{v)  +(P)  =  -  Ap  (v  -  V  ) 

^g  g  g  X  ^g  g'  g  y  '  g  y  g  p 

(C-ll) 

u  (P  )  +  V  (P  )  -  a^  u  (p  )  -  a^  V  (p  )  =  Ap  B 

g  g  X  g'  gV  g  ^g  X  g'^g  y  ^p 

(C-12) 

u  (u  )  +  V  (u  )  =  A(u  -  u  ) 
p'  p'x  p  p  y  g  P 

(C-i3) 

u  (v  )  +  V  (v  )  =  A{v  -  V  ) 

P  P  X  p  p  y  g  P 

(C-14) 

u  (h  )  +  V  (h  )  =  -  IaC  (T  -  T  ) 
p'  p'x  p  p  y  3  p  g 

(C-15) 

(u  )  dx  +  (u  )  dy  =  du 

'  g  X  g  y  ^  g 

(C-16) 

(Vg)^dx  +  (Vg)ydy  =  dVg 

(C.17) 

=  ‘^Pg 

{C-18) 

(P  )  dx  +  (P  )  dy  =  dP 

'  g  X  '  g  y  ^  g 

(C-19) 

(u  )  dx  +  (u  )  dy  =  du 
'p'x  '  p  y  '  P  • 

(C-20) 

(v  )  dx  +  (v  )  dy  =  dv 
'  p'x  '  p'y  '  P 

(C-21) 

(h  )  dx  +  (h  )  dy  =  dh 
'  p'x  p  y  P 

(C-22) 
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where 


A  = 


f 

g 

m  r  2 
P  P 


B 


=  (V  -1)  r  (u  -  a  )^  +  (v  -  V  )^  +  4  C  (T  -  T  ) 

’e  L  '  g  p'  '  e  p'  j  '  p  b'  j 


g  P 


r  -  g  Pg 
^  -  r  “F7 


The  gas  and  particle  equations  of  state  are 


{C-23) 


(C-24) 

(C-25) 


P  =  p  RT  (C-26) 

g  g  g 


T  =  f(h  )  (Tabulated)  (C-27) 

P  P 


As  discussed  in  Appendix  B,  the  characteristic  directions  thennselves 
are  found  by  solving  the  determinant 


(C-28) 


where  a.,  and  b. .  are  the  coefficients  of  the  x  and  y  derivatives  of  the 
seven  dependent  variables  in  the  seven  original  system  equations, 


L.  = 
1 


a. . 
ij 


b..u  •' 

ij  y 


+  c. 
1 


0  (i  = 


7) 


(C.29) 
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When  Equations  (30)  and  (31)  are  used,  Equation  (28)  becomes 


„  dy 


dy 


■p  (u  ^  -  V  ) 

'  g  dx  g' 


p  ,(u  -  V  )  0 

g  dx  g' 


0  p  (u  ^  -  V  )  0 

^g'  g  dx  g' 


2/„  ^ 


dy 

asT 


-1 


0  .a‘'(u  ^  -  Y  )  (u  ^-v  )  0 

'  g  d5c  g'  '  g  dx  g' 


dy 


(u  ^  -  V  )  0 

'  p  dx  p' 


<"p  ^  "V 


(C.32) 


=  0 


dv 

Expanding  Equation  (32)  and  solving  for  ^  gives 


j  2/,  u  V  +  a 


2 — r 

Ug  -  a 


dy 

u  -  V 
P  I 


=  0 


(C-33) 


From  Equation  (33),  seven  characteristic  directions  were  obtained,  the 
following  four  being  distinct: 


a)  Gas  streamlines, 


j  V 

dy  _ 

"  u 


g 


(C-34) 
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b)  Gas  Mach  lines, 

j  u  V  + 

dy  g  g  - 


-1 


•7 — r 

u  -  a 
g 


c)  Particle  streamlines, 


,  V 

dy  ^  ^ 
3x  li 

P 


(C-35) 


{C-36) 


To  determine  the  compatibility  equations  that  must  be  valid 
along  the  characteristic  curves,  the  determinant  of  the  coefficients  of 
the  14  equations  of  the  complete  system  was  modified  by  replacing  one 
column  at  a  time  by  the  nonhomogeneous  terms,  c^.  The  system  of 
Equations  (C-9)  through  (C-22)  in  matrix  form  is  shown  in  Figure  C-1. 
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All  blank  spaces  represent  zeros.  By  replacing  the  columns  one  at  a  time  with 
the  nonhomogeneous  terms,  seven  compatibility  equations  were  obtained,  each 
valid  along  a  particular  characteristic  curve.  Along  gas  streamlines: 


j  V 

dy  ^  _g 
3x  u 

g 


p  fu  du  +  V  dv  +  dP  =  -Ap  r  (u  -u  )  dx  +  (v  -V  )  dy  1 

L  g  g  g  gJ  g  ^pL  g  P  g  P 


dP  dp  Ap  Bdx 

g  _Y_  ^g  _  ^P 

F“  8  ^ - VH - 

g  '^g  g  g 


Along  gas  Mach  lines, 


,  u  V  +  a  vM^-1 
dy  _  g  g  - 

a^-  — 2 - 


Ug  -a 


(C-37) 

(C-38) 

(C-39) 


(C-40) 


(ugdy  -  Vjdx)[  AppBdx  -  UjdPj  ]  t  {  APp  [  (  »g  -  »p)  dy 


-  (''g  -  'pI  +  pg  [V"g  •  V^g'^y  ^  (V’"''g'^l  *  “V’'} '  ° 


Along  particle  streamlines, 


u  du 
P  P 


V  dv 

P  P 


u  dh 
P  P 


dy  _ 

3x  u 

p 

A(u  -  u  )dx 
g  P 

A(Vg  -  Vp)  dy 

-  4  AC  (T  -  T  )  dx 
3  P  g 


{C-42) 

(C-43) 

(C-44) 

(C-45) 
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These  seven  compatibility  equations  were  found  when  solving  for  the 
following  partial  derivatives: 


(u  )  :  (C-38)  on  gas  streamlines  and  (C-41)  on  gas  Mach  lines 

g  X 

(u^)^:  (C-38)  on  gas  streamlines  and  (C-41)  on  gas  Mach  lines 

(v  )  :  (C-38)  on  gas  streamlines  and  (C-41)  on  gas  Mach  lines 
g  ^ 

(v  )  :  (C-38)  on  gas  streamlines  and  (C-41)  on  gas  Mach  lines 

g  y 

(p  )  :  (C-39)  on  gas  streamlines  and  (C-41)  on  gas  Mach  lines 
g  ^ 

(Pg)^:  (C-39)  on  gas  streamlines  and  (C-41)  on  gas  Mach  lines 


(P  )  :  (C-41)  on  gas  Mach  lines 
g  ^ 


(Pg)^:  (C-41)  on  gas  Mach  lines 


(Up)^:  (C-43)  on  particle  streamlines 


(Up)y:  (C-43)  on  particle  streamlines 


(Vp)^:  (C -44)  on  particle  streamlines 


(Vp)^:  (C-44)  on  particle  streamlines 


(hp)^:  (C-45)  on  particle  streamlines 


(hp)^:  (C-45)  on  particle  streamlines 


Report  No.  0162-01TN-16 
Appendix  C 


When  the  flow  is  supersonic,  (M>  1),  the  system  is  totally  hyperbolic. 
The  original  system  of  eight  partial  differential  equations  can  be  replaced  by 
the  system  of  seven  total  differential  equations  valid  along  the  characteristics 
and  by  the  definition  of  the  particle  stream  function 


In  the  subsonic  flow  regime  where  M  <  1,  the  gas  Mach  lines  are  imag¬ 
inary.  However,  the  gas  and  particle  streamlines  are  still  real  and  can  be 
used  to  determine  the  gas  and  particle  properties  if  the  gas  velocity  components 

u  and  V  can  be  determined  by  some  other  procedure.  A  simplification  of  the 
S  S 

original  system  is  thus  obtained  even  in  the  subsonic  and  transonic  flow  regimes. 
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GENERAL  ONE -DIMENSIONAL  FLOW  OF  A  GAS-PARTICLE  MIXTURE 

In  the  application  of  the  nnethod  of  characteristics  to  the  supersonic  portion 
of  a  nozzle,  a  starting  line  is  required  across  the  flow  field  in  the  supersonic  re¬ 
gion.  To  determine  the  starting  line,  the  solution  for  the  entire  flow  field  in  the 
subsonic  and  transonic  portions  of  the  nozzle  must  first  be  derived.  Since  all  the 
characteristic  curves  do  not  exist  in  the  regions  with  Mach  numbers  less  than  1, 
a  different  approach  is  necessary  in  these  regions. 

The  only  feasible  approach  is  a  general  one -dimensional  analysis  in  the  sub¬ 
sonic  portion  of  the  nozzle  and  a  constant-lag  one -dimensional  analysis  in  the 
transonic  portion  of  the  nozzle.  The  general  one -dimensional  analysis  is  discussed 
in  this  appendix,  and  the  constant-lag  analysis  is  discussed  in  Appendix  E.  The 
assumptions  that  pertain  to  the  following  discussion  are  the  same  as  those  made 
for  the  axisymmetric  case  with  the  additional  restriction  of  one -dimensional  flow. 


Mass 

Momentum 

Energy 


Figure  D-1.  General  One -Dimensional  Gas -Particle  Model 


In  the  one -dimensional  control  volume  shown  in  Figure  D-1,  a  mass 
balance  results  in 


w  =  p  u  Area  (D-1) 

g  g  g 


w  =  p  u  Area  (D-2) 

g  P  P 
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A  momentum  balance  on  V  yields  the  equation 


p  u  (u  )  +  p  u  (u  )  +  (P  )  =0 

g  g  X  ‘'p  p'  p'x  g'x 


(D-3) 


An- energy  balance  on  V  results  in  the  relationship 


fc  (T  -T  )  +  4-  ^  1  +  [{h  -h  )  +  4-  u  ^  1  =  0 

®  L  Pg  g  go  2  g  J  p  p  po'  2  p  J 


1  2 


(D-4) 


where  the  subscript  o  denotes  chamber  stagnation  conditions.  The  equation  of 
state  of  the  gas  is 

P  =  p  RT  (D-5) 

g  ^g  g 


A  momentum  balance  on  a  sp^.e  rical  particle  yields  the  relationship 

u  (u  )  =3 

p'  p'x  8  m  r  (u  -u}  u  -u 

PP'g  p  ig  pi 


(D-6) 


and  an  energy  balance  on  a  spherical  particle  results  in 


u  (h  )  =  -  - (T  -  T  ) 

P  P  X  "'p  ^p  P  ® 


where 


(D.7) 


Cj3  =  f  (Re) 


Re  = 


2r  I  u  -  u  IP 
p  I  a  p  I  j 


k  Nu 


The  equation  of  state  of  the  particle  was  expressed  in  tabular  form  by  the 
relationship 

Tp  .=  f  (hp)  (Tabulated)  (D-8) 


Page  2 


Report  No.  0162-01TN-16 
Appendix  D 


When  the  system  energy  equation  (Equation  D-4)  is  solved  for  the  gas  temperature 
T  in  terms  of  the  remaining  flow  variables,  the  result  is 

O 


P 

g 


-  h  )  + 
PO 


(D-9) 


By  means  of  the  gas  and  particle  continuity  equations,  the  gas  equation  of  state, 

u  g  Area 

the  system  energy  equation  given  in  Equation  D-9,  and  the  multiplier  — a— I*- - 

w 

the  system  momentum  balance  can  be  rewritten  as  g 


W 

(M^  -1)  (u  )  +  — 

g  X 

g 


7^—  u  )(u  ) 
C  p  p  X 
Pg  ^ 


RT. 


Area 


^  (Area) 


=  0 


where 
2 


M 


2 
u 


As  in  the  axisymmetric  case. 


f  = 


Stokes 


Nu 

ITTiry 


Cp 

24 

Nu 


Stokes 


the  following  definitions  apply: 
and 


(D-IO) 


(D-11) 


(D-12) 


Equations  D-6  and  0-7  can  then  be  solved  for 


(u  )  and  (h  ) 

'  p'x  '  p'x 


du 

u  -  u 

A  (  -£ - E.) 

u 

(D-13) 

p 

35^  = 

2 

-  -j-  AC  (  — E_ - B.)  where 

(D-14) 

P 
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A  = 


9 


m  r 
P  P 


and 


(D-15) 


C  -  ® 
C  -  j- 


Pg 

5r 


.  (D-16) 


Substituting  Equations  D-13  and  D-i4  in  Equation  D-10  results  in 


du 

g  _ 

“35?  ■ 


u 

g _ 

Area 


d(Area) 

3x 


■)  -  1^1  (H 

I  Of  '  >  tJ  '  u 


+  A 


^  “r.  j“g-"p)}  ]] 


a  =  y  RT 
g  g- 


(D-17) 


Solving  Equations  D-1  and  D-2  for  p  and  p  results  in  the  following 

g  P 

relationships 

w 

g 
P 


g  '  Ug  Area 


p  u Area 
^  P 


(D-IS) 

{D-19) 


The  gas  pressure  is  found  from  the  gas  equation  of  state 

P  =  p  R  T 
g  g 

The  shape  of  the  nozzle  determines  the  function 
Area  =  £  (x) 


(D-20) 


(D-21) 


the  cross-sectional  area  of  the  nozzle. 
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Thus,  the  general  one -dimensional  flow  of  a  gas-particle  system  is 
represented  by  the  following  system  of  equations. 


du  *  u  -u 
P  =  A  g  .  p 

dx  u_ 


dh 

asT 


.  =4  (!^) 
u 


M^-1  I  t 


d  (Area) 
area  332 


iT  =  y  R 

g  g 


a 

Area  =  f  (X) 


Z  - 2 

m  r 
P  P 


f  =  f  (Re) 
g  =  g  (R«) 


g 


Pg  “  u  Area 
g 

^P 

**0  “  u_  Area 


(Given) 


(Given) 


(Given) 


(D-22) 


(D-23) 


(^)  fz  [I  (v,-‘)  AC  (r^.T^) 

+  A  <1  Y  (u  -u  )^  +  u  (u  -u  )\  1 
\  ^g  g  P  P  g  P  J  -I 


(0-24) 


(D.25) 

(D-26) 

(D-27) 


(D..28) 


(D-29t 
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Tp  =  f(hp)  (Tabulated) 


(D-30) 

(D-31) 


P  =  p  R  T 
g  ^g  i 


V  =  f  (T  )  and  c  =  f  (T  J  ( Tabulate d.) 

^g  g  Pg  g 


H  T 

g  o  g 


(D-32) 

(D-33) 

(D-34) 


In  Equations  D-22  through  D-34,  the  quantities 

Pr,  nip.  Tp,  Wg,  Wp,  hp^,  R,  a  are  all  given  constant  parameters . 


Equations  D>22,  D-23,  and  D-24  can  be  solved  numerically  in  the 

subsonic  regime  for  u  ,  h  ,  and  u  by  any  of  several  standard  techniques, 

P  P  g 

such  as  the  Adams'  or  the  Range -Kutta.  method,  provided  that  starting 

values  are  assumed  at  some  location  where  the  flow  velo  '^ities  are  small 

and  gas -particle  equilibrium  can  be  assumed.  The  remaining  flow  variables, 

p  ,  p  T  ,  T  .  and  P  ,  then  can  be  calculated  from  Equations  D-28  through 
g  p  I  .p  g  g 

D-32,  respectively.  Thus,  the  entire  one -dimensional  subsonic  flow  field  is 


determined  up  to  some  value  of  M<  1.  Near  M 

1 


1,  the  n’imerical  solution 


begins  to  diverge  as  a  resvilt  of  the 


M^-1 


term  in  Equation  D-24. 


Because  of  this  divergence,  it  becomes  necessary  to  consider  a  constant 
lag  one -dimensional  gas -particle  system  to  obtain  solutions  in  the  transonic 
regions. 
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To  obtain  a  transonic  solution  for  a  one -dimensional  gas -particle  mixture, 
it  was  necessary  to  assume  that  the  particle  kinetic  and  thermal  lags  are  con¬ 
stant  throughout  the  transonic  regime.  The  particle  kinetic  and  thermal  lags  were 
defined: 

(E-1) 
(E-2) 

It  was  also  assumed  that  the  particles  do  not  change  phase  in  the  transonic 
portion  of  the  nozzle.  The  particle  equation  of  state  was  then  expressed  as 

‘‘■p  -  V’  “  'pi  '’■p  • 

where  in  the  chamber, 

go  po 

The  temperature  of  the  gas,  was  found  to  be 

(E-4) 


(E-5) 


Substituting  Equations  (E-l),  (E-2),  and  (E-3)  in  Equation  (E-4)  results 
in  the  following  expression  for  T^: 

T  -T  !IiL 

g  ^go  7 


u 

K  =  ^  0  <  K  ^  1 
g 

T  -T 

L  =  0  L  ^  1 

go  ‘  g 
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The  particle  energy  balance  was  fovind  to  be 


A  C 


Substituting  Equations  (E-1),  (E-2),  {E-3),  and  (E-5)  in  Equation  (E-12) 
gives 


(E-12) 


du  , 
g  -  1 

asr  ■■  T 


A  c 


(l-L) 
c  ,  K  L 
14 


Equating  Equations  (E-11)  and  (E-13)  gives 


(E-13) 


L  = 


1 


1  + 


{E-14) 


Equation  (E-14)  determines  the  relationship  between  the  particle  thermal 
and  kinetic  lags  in  a  constant -lag  situation. 


Eqtiation  (E-5)  was  solved  for  the  ratio  T  /T  to  give 

9  9 


g 


Yg  =  ‘  *  Wg  -1)  B 


Equation  (E-7)  was  considered  at  the  geometric  throat  where 


and 

Therefore 


d(Area)  _  d  (Area)  dx  _  ^ 

(Area)  “  dx  (Area) 


du  du 


dx 


>  0 


=  ^  22 
dx  u„ 
g  g 

RT  j  -  (ij  =  0  and  =  1  at  the  throat. 


(E.i5) 

(E-16) 


(E-17) 
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Therefore,  from  Equation  (E-15)  at  the  throat,  where  conditions  are 
denoted  by  *, 


T*  =  ^  T 

«  (Vg  +  1) 


(E-18) 


Equation  (E-7)  was  integrated  in  the  following  manner: 

-  -1)  ^  (Equation  E-7) 


_  ?  ,  Eu 

g 


1  g 


T  =  T  -  -  - 

g  go  Z  c 


M  ^ 


Eu 


Pg 

2 


_ g _ 

u  2 

V  R(T_-  D) 


'g  go  2  c 


pg 


d(Area)  _ 
( Area)  ~ 


Eu  du 

_ g  g 

j  u  2 

y  R(T  -  i 
'g  '  go  2 


pg 


du. 


u 


D) 


In  A  =  -  ^ - jy  jj  In  Tg  -  In  u^  +  Constant 


(E-19) 


Evaluating  the  above  expression  between  a  general  point  and  the  throat 
*  — 

where  and  M  =  1  gives 


(Area) 
(Area)  * 


(l  +  ^ 

\  2 

(E-20) 
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The  gas  continuity  equation,  the  perfect  gas  law,  and  the  system 
momentum  equation  were  combined  to  determine  the  following  pressure 
relationship: 


'i/  =  p  u  Area 

g  ^g  g 


(E.21) 


^^g  I  ^'^g  ,  d(Area) 


0 


P 

g 


RT 


dP 


dT 


g 


d(Area) 

(Area) 


-1) 


du 

-Jt 

"g 


Combining  the  above  gives 


(E-22) 


(E.23) 


dP  2  du  dT 

.m"  ^  -  ^  -0 

g  g  g 


In  the  same  manner  as  Equation  (£-7)  was  integrated,  Eq\iation  (E-24) 
was  solved  and  evaluated  between  the  general  condition  and  the  stagna¬ 
tion  chamber  condition  where  P  =  P  and  T  =  T  to  give 

g  go  g  go  ® 


(E-24) 


(E-25) 
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Substituting  Equations  (E-15)  and  (E-25)  into  the  perfect  gas  law  gives 


(E-26) 


Combining  the  gas  and  particle  continuity  equations  gives 


P 


P 


1 

K 


P 


g 


The  solution  of  a  one -dimensional  gas -particle  system  under 
constant -lag  conditions  is  thus  given  by 


(E-27) 


(Area)  1  /  2 

■  ra  V  + 1 

'  •  or 


2(Vg  -1) 


{E-28) 


g 


(E-29) 


(E-30) 


(E-31) 
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M  ^  =  E  M  ^ 


2  xx2  2 

u  =  M  a 

g 


a  =  Y  RT 
^g  g 


P  g 


T  ={1-L)T  +LT 
P  go  g 


i  w 
^  P 

P  ^  —  '8 

g 


1 


1  +  3 


t 


'g 


D  = 


1  +  (Vg  '  1)  X 


<v 

K 

8 _ 


w  c  , 

it-E-  -ELl 

C  j 

g  Pg  J 


E  = 


(‘-■^1  *  ■']"  K-‘l 


pl 

pg 


L  D 


} 


(E-32) 

(E-33) 

(E-34) 

(E.35) 

(E-36) 

(E-37) 

(E-38) 

(E-39) 

(E-40) 


(E.41) 


Equations  E-28  through  E-31  are  seen  to  be  the  one -dimensional  isentropic 

gas -dynamic  relationships  withY  sind  M  replaced  by  the  modified  parameters 

__  g 

and  M,  respectively.  A  gas -particle  mixture  in  a  constant-lag  nozzle  with  no 
phase  change  can  thus  be  considered  as  a  perfect  gas  with  appropriately  modified 
specific  heat  ratio  and  Mach  number. 
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From  Equation  (E-11)  it  is  seen  that,  for  constant  lag  nozzles, 

A  (1  -  K)  ^  ^ 

— s.  =  - ! — , — L  =  Constant 

QX  ■fm2. 


Thus,  the  assumption  of  constant  lag  imposes  the  restriction  of  constant 
axial  velocity  gradient.  Such  a  condition  is  approximated  in  the  throat  region  of 
most  nozzles.  The  solution  for  constant-lag  flow,  therefore,  shows  that  in  the 
throats  of  most  nozzles  the  gas -particle  mixture  may  be  treated  as  a  perfect  gas 
with  the  appropriately  modified  specific  heat  ratio  and  Mach  number. 
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TRANSONIC  FLOW  APPROXIMATION 


By  using  the  analysis  of  a  one -dimensional  gas -particle  mixture  under 
constant-lag  conditions,  it  was  found  that  the  gas -particle  mixture  could  be 
treated  as  a  perfect  gas  with  appropriately  modified  specific  heat  ratio  and 
Mach  number.  The  solutions  for  transonic  flow  can  therefore  be  obtained  in  the 
same  manner  as  for  a  perfect  gas,  since  the  two-dimensional  gas -particle  flow 
in  the  throat  region  can  be  treated  approximately  as  a  one-dimensional  flow,  and 


the  actual  specific  heat  ratio  and  Mach  number  of  the  two-dimensional  flow  can  be 
approximated  by  the  modified  parameters  and  ^  based  on  the  one -dimensional 
approximations.  Several  methods  have  been  proposed  for  analyzing  transonic 
perfect-gas  flow.  Methods  proposed  by  Sauer,  Oswatitsch  and  Rothstein,  and 
Kliegel  (References  F-1,  F-2,  and  F-3,  respectively)  appear  to  be  practicable, 
Sauer's  method  will  be  used  until  a  better  solution  is  available. 


Figure  F-1  illustrates  the  geometry  of  the  model  considered  by  Sauer, 


Figure  F-1.  Transonic -Flow  Model 
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The  assumptions  made  for  the  model  are  as  follows: 

1.  The  gas  is  a  perfect  gas. 

2.  The  flow  is  irrotational. 

3.  The  flow  is  two-dimensional  (6  =  0)  or  axisymmetric  (  6  =  1). 

4.  The  specific  heats  are  constant. 

By  combining  the  continuity,  momentum,  and  speed-of-sound  equations  for 
such  a  flow  field,  the  following  equation  was  obtained: 

~/2~2.  ~,2~2,  -'-~~..2v  - 

u  (a  -u  )  +  v  (a  -V  )-2uvu  +6a  —  =0  (F-i) 

X  y  y  y  ' 

where  u*  and  “v  are  the  x  and  y  components  of  velocity,  and  a  is  the  speed 
of  sound.  The  energy  equation  is 


1/^2  2 , 

2-(u  +v  )  +  h  =  constant 


{F-2) 


~  2  »«»  2 

u  +  V 


+  a  =  constant 


{F-3) 


^2  ^2  2 

Defining  o  as  the  speed  of  sound  where  (u  +  v  )  =  a  ,  Equation  (F-3) 
becomes 


a  = 


2  I  Yo  +  M  1 


'g  ■  I  -  *2  »g 
—2 -  *  \  — 2" 


2  M  2 

u  +  V 


(F-4) 


The  dimensionless  velocity  components  are  defined  by 


Tr  =  i,  and  V  = 


(F-5) 


Page  2 


Report  No.  0162-01TN-16 
Appendix  F 


In  terms  of  U  and  V,  (F-1)  becomes 


U 


+  1 


uvu 

y 


I 


Y  jT  "  ^  2  2 

-  (TT^  +  =  0 

*  ‘ 


(F-6) 


By  limiting  the  investigation  to  the  vicinity  of  the  sonic  condition  where  U  ^  1  and 
V  ^  0,  TJ  and  V  can  be  expressed  in  terms  of  perturbation  velocities  as 


U  =  1  +  u  and  V  =  v 


(F-7) 


where  u  and  v  are  much  smaller  than  1.  In  terms  of  u  and  v.  Equation 
(F-6)  becomes 


u 

X 


Y  *  i 

•  nr 


'  Y  +  1 

g 


+  -  u  (1  +  u)  v 

Y  +  1  ^ 


Y  +  1 


-  2 


Y  +  1 


r  ^ 

.'.2  r8'‘i 

1 

pi 

1 

L  Y  +  1 
g 

•  Y  +  1 ' 

+  V 


1 


J 


+  V 


'] 


=  0 


y 

(F-8) 
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As  X  and  y  approach  zero,  u  and  v  approach  zero.  However,  should  not 
approach  zero.  Equation  (F-8)  shows  that  as  x  and  y  approach  zero,  ^  and 
approach  zero.  From  the  symmetry  about  the  x  axis,  u^  approaches  zero 
at  y  =  0.  Substituting  these  limiting  values  into  Equation  (F-8)  and  neglecting 
terms  of  higher  order  than  the  first  in  quantities  which  approach  zero  as  x  and 
y  approach  zero.  Equation  (F-8)  becomes 


(Vg  +  1 )  u  u^ 


{F-9) 


For  irrotational  flow, 

Qu  _  8v  (F-10) 

ay  ~  ax 

Hence,  a  velocity  potential  <j)  (x,y)  can  be  defined  such  that 

u  =  d  and  v  =  ?(  (F^ll) 

The  velocity  potential  can  be  defined  as  a  power  series  in  y,  containing 
only  even  powers  of  y  because  u  is  symmetrical  about  the  x  eixis. 

d  (x,y)  =  f^  (x)  +  f^  (x)  y^  +  f^  (x)  y^  +  ...  (F-12) 

u  (x,y)  =9»^  =  f^(x)  +  f2'(x)y^  +  f4'(x)y‘‘+  ...  (F-13) 

V  (x,y)  =  =  Zf^  (x)  y  +,  4f^  (x)  y^  +  .  .  .  (F-14) 


where  the  primes  denote  ordinary  differentiation  with  respect  to  x. 
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Substituting  the  power  series  for  u  and  v  into  Equation  (F-9)i  equating 
the  sum  of  like  powers  of  y  to  zero,  and  solving  for  f^  (x)  and  f^  (x)  gives 


(X) 


<v 

2  (3  +6) 


(F-15) 


(F-16) 


When  y  =  0i  u  (x,  0)  =  f^  {x),  where  u  (x,  0)  is  the  axial  perturbation 
velocity  distribution  along  the  x  axis.  If  u  (x,  0)  is  known,  f^  (x)  and  f^ 
(x)  can  be  determined,  and  the  flow  field  would  be  established.  For  a 
linear  axial  perturbation  velocity  distribution. 


u  (x,  0)  =  a  X 


(F-17) 


+  1)  X 


f-  (x)  =  - 2 - 

2  '  2  (1  +5) 


^4  (x)  = 


..2  3 

(v  +  1)  a 

8(1+5)  (3  +5  ) 


(F.^8) 


(F-i9) 


,  .  +  2 

u(x,y)=  ax  +  — 2 -  y  +  ... 

2  {1  +  5) 


(F-20) 


(y  +  1)  a  X  (y  +  1)  a  ^ 

v(x,y)=;  — 2 -  y+  _S -  y  + 

(1  +5  )  2(1  +5  )  (3  +6  ) 


.  .  (F.  2r. 


—  2  2  2 

The  critical  curve  where  M  =  1  and  (  u  +  v  )  =  a*  can  be 
determined  as  follows: 


«,2,«,2  —  Z  *2,-2  *2 

u  +  V  =  a*  =  U  a*  +  V  a* 


(1  +  u)  ^  +  V  ^  =  1 
u  =  0 


iF-ZZ) 
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Therefore,  the  critical  curve  is  found  from  Equation  (F-20)  to  be 

(Vg  +  1)  “  2 

— - -  y 

^  "  ■  2  (1  +  6  )  (F-23) 


To  locate  the  coordinate  system  in  the  nozzle,  let  y  at  x  =  0  be  equal 

to  y  ,  and  solve  for  6  from  Eqxiation  {F-23); 
s 


+  1)  a 

— -  y 

2(3+5)  ® 


{F-24) 


e  is  the  distance  from  the  geometrical  throat  downstream  to  the  gas 
dynamic  throat  where  M  =  1. 

Next,  by  using  Figure  F-2,  determine  the  curvature  of  the  nozzle 

wall  at  the  narrowest  cross  section.  Curvature  is  defined  as  the  change 

of  direction  of  the  tangent  to  a  curve  per  unit  distance  along  the  curve, 

and  the  radius  of  curvature  is  denoted  by  p  . 

0  ® 


Figure  F-2.  Radius -of-Curvature  Model 
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The  flow  field  near  the  sonic  line  in  the  throat  of  a  nozzle  is  thus 
determined  by  the  following  relationships: 

{\g  +  1)  ^  2 

u  (x,y)  =  a  X  +  - = -  y 

2(1+6) 


(F-28) 


V  (x,  y) 


(y  +  1)  a  ^  X  (y  +  1)  ^  a  ^ 

— S -  y  +  _£ - 


(1  +6) 


2  (1  +  6  )  (3  +  6  ) 


(F-29) 


2  (3  +  6  ) 


(\g  +  1)  (1  +*>  )  Yg 


{F-30) 


/  (1  +t>  ) 

Ps  ^s 


(F-31) 


where  y  and  p  are  geometric  properties  of  the  nozzle, 
s  s 


By  using  the  modified  specific  heat  ratio  and  Mach  number  defined 
for  gas -particle  flow  under  constant  lag  with  linear  velocity  gradient  and 
no  phase  change,  the  above  equations  can  be  used  to  determine  the  two- 
dimensional  or  axisymmetric  flow  of  a  gas -particle  system  in  transonic 
flow, 

u  (x,  y)  =  a*  (1  +  u) 


(F-32) 


V  (x,y)  =  a  V 


(F-33) 


2^  R  T, 


(F-34) 
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u. 


V 

g 


_1 

(E)  u 


(F-35) 

(F-36) 


where  E  and 


are  defined  in  Appendix  of  this  report. 
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DERIVATION  OF  THE  FINITE -DIFFERENCE  EQUATIONS 


Seven  total  differential  equations  were  derived  that  are  valid  along  the  four 
characteristic  curves.  Because  of  the  nonlinear  nature  of  these  equations, 
numerical  techniques  had  to  be  used  to  obtain  a  solution.  The  seven  compatibility 
equations  were  put  into  finite-difference  form  and  are  presented  in  this  appendix. 
The  remaining  equations  necessary  to  carry  out  the  numerical  solution  by  use  of 
a  high-speed  digital  computer  were  also  put  into  a  form  suitable  for  computation. 


Three  equations  for  particle  properties  were  derived  for  a  constant  par¬ 
ticle  size.  To  account  for  the  possibility  of  more  than  one  particle  size  or 
species,  six  discrete  particles  were  allowed  for  in  the  numerical  solution.  This 
was  accomplished  by  considering  the  streamline  of  each  additional  particle  as  an 
additional  characteristic  curve  and  by  applying  the  particle  compatibility  equations 
independently  for  each  particle  along  its  streamline. 


The  equations  were  put  into  a  form  suitable  for  calculation  of  an  interior 
point  in  the  flow  field.  For  special  points,  such  as  the  nozzle  axis,  the  nozzle 
wall,  the  limiting  particle  streamline,  and  a  free  boundary,  the  same  equations 
apply;  but  the  procedure  for  obtaining  a  solution  must  be  modified.  The  character¬ 
istic  net  for  such  an  interior  point  is  bounded  by  the  initial  data  line  and  the  two 
Mach  lines,  as  illustrated  in  Figure  G-1. 
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Figure  G-i.  Interior- Point  Characteristic  Net 

In  Figure  G-1,  functions  evaluated  along  gas  Mach  lines  are  designated  by  sub¬ 
scripts  I  and  II,  corresponding  to  right  and  left  running  Mach  lines,  respectively; 
and  functions  evaluated  along  gas  and  particle  streamlines  are  designated  by  G 
(gas)  and  s  (solid),  respectively. 

Along  right  and  left  running  Mach  lines,  the  characteristic  and  compatibility 
equations  are: 
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To  integrate  these  equations,  all  properties  were  assumed  to  be  constant  at  the 
average  value  of  the  property  along  the  line  of  integration.  Equations  G-l  and 
G-2  then  became 


I  .  u  V  ±  a  - 1 

=  — rr“ 


(G-3) 


+  QAu^  +  RAv  =S 

g  g  g 


(G-4) 


where  the  bar  denotes  average  values,  and 


a  =  R  T_, 

g  g 


_  Z  _  2, 

-  2  +  ''a  ) 

M  *  — S — — 
a  ^ 


M  = 


a  A  y  -  Ug  K 


0  =  5^  X  V  A  X 
g  g 


R  =  -  a  ^  p  u  Ax 
^g  g 


S  ={(y 


P  V  - 

^8  8  t2t 


-2 


y  d 


i=l 

1  i=l 

K  = 

K 

Ay  -  Axj 

A  = 

9 

7 

m  r 

P  P 
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B  =  (Vg  -  1)  [  (Ug  -  +  (Vg  -v/  +  I  C  (Tp  -  Tg)  ] 


=  f(Tg) 


C  =-i_  -Pi- 

T  Pr 


'Pg  -  f  (Tg) 


Tabulated 


f  =  f  (Re) 


g  =  g  (Re) 


2r  p 
P  g 


(u  -  U  )^  +  (v  -  V  )  ^ 

g  p  g  p 


Pr  =  Constant 


U  =  T 


where  u  and  a  are  constants, 

M3 


P  =  p  R  T 
g  '^g  g 


(G-5) 


Since  Equation  (G-4)  is  valid  along  both  Mach  lines,  two  equations 

relating  AP  ,  A  u  ,  and  Av  were  obtained: 
g  g  g 


t 

gl 

"  “I 

Au  ^ 
gl 

AV  ^  =:  S. 

gl  I 

(G-6) 

gll 

^  ^gll 

+  Rj, 

A  V  . =  S  -T 
gil  il 

(G-7) 
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Along  the  gas  streamlines,  the  characteristic  and  compatibility 
equations  are 


'dy  _  ^ 
ax  ^  u 

1  /  G  g 


{G-8) 


p  r  u  du  +  V  dv  I 

^gL  g  g  g  g  J 


I  [ 


6 

"E  {  <'^g-''p)dy  ]  }  (g-9) 

i=l 

AppBdx  (G-10) 

“FH - 

g  g  J' 


In  integrated  form,  these  equations  are 


Ij^ 

1 

V 

=  -1- 
G  u„ 

(G-ll) 

g 

A  P 

gG  G 

(G.12) 

^PgG  = 

-h  "Tc 

(G-13) 

U  '^g 

V  =  p  V 

^g  g 

6 

■  Z  { [ <"g"  V  ^^g  ■  V  ^ ^  ]  } 

i=l 


V  F  u 

'g  g  g 


) 


A  X 
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Along  the  particle  streamlines,  the  characteristic  and  compatibility 
equations  for  each  particle  are 


( 


u 

p 


(G-14) 


u  du  =  A  (u  -  u  )  dx 
P  P  g  P 


V  dv  =  A  (v  -  V  )  dy 
P  P  g  P  ' 


(G-15) 

(G-16) 


u  dh  =  -  4-  AC  (T  -  T  )  d  X 
P  P  ^  P  g 


(G-17) 


In  integrated  form,  these  equations  are 

^ 

S  u 

P 


'  u  -  u_ 

A  =  A  M- 


u 


(G-18) 


(G-19) 

(G-20) 

(G-21) 
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The  particle -stream  function  equation,  when  solved  for  p  ,  becomes 

P 


Pn  = 


P  “  c(Updy  -  Vpdx) 


In  integrated  form,  this  equation  becomes 

y-'l'pl  ■  "pi  ■  ^2  " 


"p4  - 


1  '^Pl 

TET'VTjr 


£1 


(G-22) 


Li  =(yj 


u  {y4-yi) 


^2  =  ®2 


(y4  ■  ^2^  •  %  ^*4  ■  *2)> 


'^p4  “^p2  +  ^2  ^Pp2  ■'■Pp4) 


where,  for  two-dimensional  flow. 


»!  =  Oj  =  I 


and,  for  axisymmetric  flow. 


1 


■l  ■  7  '^1  *  ''4' 


"z  =  7  <''2  *  ’'4> 
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Equations  G-6,  G-7,  and  G-12,  relating  the  gas  pressure  and  velocity 
components,  can  be  solved  for  these  properties  at  the  intersection  of  the  Mach 
lines  to  give 


(Y-MjX)  (Rjj-Mjj  Vq)  -  (Z-MjjX)  (Rj-Mj  7^ 

Ug4  =  _  ^ 

(Qj  -Mj  TT^  (Rjj  -MjjVq)  -  (Qjj  -  (Rj  - 

(Qj  -  MjTJ^)  (Z-Mjj  X)  -  (Qjj  -  (Y  ~  MjX) 

V  "  (G-24) 

(Qj  -  Mj^q)  (Rjj  -  ’  ^I^G^ 


P.=X-U'^u.-V„v.  +  P- 
g4  G  g4  G  g4  g3 


(G.25) 


X  =  D_  +  V_  u  -  +  T„  V  , 
G  G  g3  G  g3 


Y  =  s,+M,  (Pg,  -Pg3)  +  Q:u^,  +  Rivg, 


Z  =  Sjj  +  Mjj  (P^2  -  Pg3>  +  °II  "gZ  +  Pu  ''gZ 


The  gas  density  at  point  4  from  Equation  G-13  is  given  by 


g4  ^g3 


+ 


y  R  T 
''g  g 


+  T, 


where  T_  is  defined  above. 

(j 


(G-26) 
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Equations  (G-19))  (G-20),  and  (G-21)  were  solved  for  the  particle  velocity 


components  and  enthalpy  at  point  4  to  give 

(G-27) 

(G-28) 

V  =  ^5  + 

(G-29) 

The  particle  density  at  point  4, 


p4’ 


is  given  by  Equation  (G-22). 


Thus,  the  original  system  of  characteristic  total -differential  equations 
has  been  put  into  finite-difference  form,  and  expressions  for  all  the  flow 
properties  at  point  4  have  been  obtained  in  a  form  suitable  for  computation. 
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NUMERICAL-SOLUTION  PROCEDURE  FOR  COMPUTER  PROGRAM 

I.  GENERAL 

The  equations  governing  a  gas -particle  system  have  been  derived,  the 
characteristic  equations  have  been  determined,  and  the  characteristic  equations 
put  into  finite-difference  form  in  the  preceding  appendixes.  In  this  appendix,  the 
use  of  these  equations  is  demonstrated  for  determining  the  flow  properties  at  an 
interior  point  in  the  flow  field:  any  general  point  that  is  not  on  a  solid  or  free 
boundary,  on  the  nozzle  axis,  or  in  the  region  between  the  limiting  particle  steam¬ 
line  and  the  nozzle  wall.  For  these,  slight  modifications  to  the  procedure  are 
necessary. 

This  procedure  has  been  programed  for  the  IBM  7090  computer.  The  com¬ 
puter  program  follows  the  solution  procedure  outlined  in  this  section. 

The  units  of  the  physical  parameters  appearing  in  the  solution  are  shown  in 
Table  H-1.  The  given  data  and  the  solution  are  both  expressed  in  these  units. 

Two  types  of  data  are  needed  for  this  solution:  general  and  specific.  The 
general  data  (Table  H-2)  pertains  to  the  entire  nozzle  calculation,  and  the 
specific  data  (Table  H-3)  pertain  only  to  the  flow  properties  at  the  point  under 
consideration.  The  properties  found  at  point  4  consist  of: 

1 .  Location  of  point  4  (x  and  y) 

2.  Gas  properties;  the  same  as  the  initial  data 

3.  Particle  properties;  the  same  as  the  initial  data. 

The  numerical  solution  is  based  on  the  characteristic  net  shown  in  Figure 

H-1. 

II.  SOLUTION  PROCEDURE 

The  detailed  solution  procedure  for  the  location  of  point  4  and  the  calcula¬ 
tion  of  the  flow  properties  at  that  point  is  in  the  following  sequence: 
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Re 


TABLE  H-1 

UNITS  OF  PHYSICAL  PARAMETERS 

specific  heat  at  constant  pressure 
specific  heat  ratio  of  the  gas 
ratio  of  to  for  Stokes'  flow 

ratio  of  N^  to  N^  for  Stokes'  flow 
particle  enthalpy 

particle  density  per  unit  volume  of  particle 
gas  static  pressure 
Prandtl  number 

particle  radius,  based  on  spherical  particles 

gas  constant 

Reynolds  number 

gas  static  temperature 

particle  static  temperature 

axial  gas  velocity  component 

axial  particle  velocity  component 

normal  or  radial  gas  velocity  component 

normal  or  radial  particle  velocity  component 

axial  coordinate 

normal  coordinate 

component  of  gas  viscosity  expression 
gas  viscosity 

coefficient  of  gas  viscosity  expression 
gas  density 

particle  density  per  unit  volume  of  gas 


Btu/lbm-  °R 
dimensionless 

dimensionless 

dimensionless 

Btu/lbm 

Ibm/ft^ 

Ibf /in.  ^ 
dimensionless 
micron 
Btu/lbm- “R 
dimensionless 

*R 

ft /sec 

ft/sec 

ft/sec 

ft/sec 

ft 

ft 

dimensionless 
Ibm/ft-sec 
lbm/ft-sec-(  ®R)° 
Ibm/ft^ 

Ibm/ft^ 
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TABLE  H-2 

GENERAL  INITIAL  DATA 

1.  Two-dimensional  or  axisymmetric  flow  specification  (5  =  0  or&=  1) 

2.  Number  of  particles  and  the  density  m^,  and  radius  r^  of  each 

3.  Tabulated  values  of  f  as  a  function  of  Re 

4.  Tabulated  values  of  g  as  a  function  of  Re 

5.  Tabulated  values  of  y  ,  R,  C  .  U  .  and  Pr  as  a  function  of  T 

g  Pg  g  g 

6.  Tabulated  values  of  T^  as  a  function  of  h^  for  each  particle  species 

7.  The  solidfication  temperature,  T^^,  of  each  particle  species 

8.  The  iteration  accuracy  criteria,  K  and  L 
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TABLE  H-3 


SPECIFIC  INITIAL  DATA 
(given  at  both  points  1  and  2) 


1.  Location  of  points  1  and  2,  (  x  and  y). 

2.  Gas  properties  at  points  i  and  2: 


-  axial  velocity  component 

-  normal  velocity  component 

-  density 

-  pressure 

-  temperature 


3.  Particle  properties  for  each  particle  at  points  1  and  2: 


axial  velocity  component 

normal  velocity  component 

density 

enthalpy 

temperature 

stream  function 
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Particle  Streamlines 
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1.  Determine  the  location  of  points  1  and  2,  and  the  gas  and  particle 
properties  at  points  1  and  2: 


V 

u  .  ,  u  , 
pi  p2 

(x^.  y^) 

'"gl’  '^82 

pi  p2 

(X2.  Y2) 

Pgl’  Pg2 

^p 1 ’  ^p2 

^gl’  ^g2 

h  .  ,  h  - 
pi  p2 

Jgl  ’  '^82 

pi  p2 

LV’  pu 

2.  Determine  T  at  the  midpoint  of  line  1-2,  and,  from  tables  of  v 
8^^  8 
and  C  vs  T  ,  determine  c  and  v  at  this  temperature.  This  c  and  v  are 
Pg  8'  Pg  8  Pg  8 

then  considered  as  the  constant  values  of  c  and  v  for  the  entire  net. 

pg  ^8 


3.  Calculate  u  ,  v  ,  u  ,  v  ,  and  p,  at  the  midpoint  of  line  1-2, 
8  8  P  P  8 


m 


-2-(Pi  +P2) 


where  P  is  any  of  the  above  properties, 
point  of  line  1  -2. 

^^8m  = 


Calculate  the  gas  viscosity  at  the  mid- 

U  T  ° 

^o  gm 


4.  Calculate  Re  for  each  of  the  6  particles  at  the  midpoint  of  Line  1-2. 

Re  =  6.5602  X  10“^  J  r -\/(u  -u  J  +  (v  -v  I 

»^gm  1  P  V  gm  pm'  gm  pm'  j 

From  tabulated  values  of  f  and  g  vs  Re,  f  and  g  for  each  of  the  particles  may  be 
determined.  Calculate  A  and  C  of  the  particles. 
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li 


A  -  4.  1806  X  10  m  r  Z" 

P  P 

C  =  25=cl0^f 

5.  In  the  following  steps,  note  that  the  value  of  any  flow  property  P  at 
point  i  on  the  initial  data  line  1-2  can  be  found  as  follows: 


iX.  -  X,  \ 


6.  Locate  point  4.  For  the  first  overall  approximation  of  point  4,  base 
the  location  on  properties  at  points  1  and  2  only.  For  subsequent  approximations, 
base  the  location  on  the  average  of  properties  at  points  1,  2,  and  4. 


ai  = 


2  z"!- 

'^81  ""81  * 


a,  =  32.  139  y  RT  , 

2  ^g  g2 


2  ^  2 
%2  +  %2  « 


®gi 


®g2  '  ’  '''g2'"g2’ 


°gl  '  ‘  (1/Mj) 


°g2  ^ 


^11  ~  ^®gl  ■  ®1^ 


^112  ^  •®g2  °2^ 


For  the  first  overall  trial  on  point  4,  set  =  pjj  andpjj^  - 


On  subsequent  trials,  and  will  be  known  from  flow  properties  at  point  4. 
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^^14  ■  ~r  ^  ■*■^^14^  ^^24  ■  T  ^^^112  ■'■^^114^ 


m  .  =  tan  p 


^Z4  =  >^24 


^2  ”^24  ~  ^1  "^14  +  <yi  ~  ^2^ 
(^24  •  ^14^ 


^4=^2  +  ”"24  ^^4  -  ^2^ 

■  ^  u  i 

7.  Locate  point  3,  For  the  first  overall  trial  for  point  4,  =  -j 

{yLli2  +  IJLii.).  On  subsequent  trials,  will  be  known  from  flow  properties  at 

point  4.  Calculate  During  the  first  overall  trial  on  point  4,  assume 

X  =  T  (x  +  X,)  as  a  first  approximation.  During  subsequent  overall  trials  on 

3  to  1  to 

point  3,  assume  that  x^  is  the  final  value  of  X3  from  the  previous  overall  trial  on 
point  4. 

Pg3  =  tan  <'^g3/'"g3)  ^^C4  = 


^^34  =  T<^^3  ^G4> 


""34  =  >^34 


Pass  a  line  through  point  4  with  the  average  slope  and  determine  a  new 

location  for  point  3  as  the  intersection  of  this  line  with  1-2.  Note  that 
”"12  ""  ■  "^2^  constant. 


1 


”"12  *2  ■  ”"34  * 


^4> 


73  =  +mj2  (X3  -  *2) 
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Repeat  the  above  steps,  obtaining  even  more  refined  values  of  and  each 
time  using  the  last  location  of  point  3  to  calculate 


Distance 


3(n  +1)  -3( 


HJIstance  f  1  -i] 


ai 


where  K  is  an  accuracy  criterion  and  n  denotes  the  nth  trial.  When  the  above 
equality  is  satisfied,  point  3  is  located. 


8.  In  all  the  equations  which  follow,  all  flow  properties  must  be 
expressed  as  average  values.  This  is  accomplished  by  replacing  each  property, 

P,  by  the  relationship 

■p  =  4-  (p,  +  p.^) 

where  P.  is  the  value  of  the  property  at  point  i  and  P^^  is  the  value  of  the  property 
at  point  4  where  F  is  evaluated  along  line  i-4.  On  the  first  over-all  trial  on  point 
4,  all  values  of  P.^  are  set  equal  to  P^.  On  subsequent  trials,  actual  values  of 
properties  at  point  4  will  be  used. 


9.  Set  two-dimensional  or  axisymmetric  flow  parameters. 


If  d 

y 

II 

0 

P 

Q 

II 

0 

SIGK 

=  0.0 

SIGK  =1.0 

SIGl 

=  1.0 

SIGl  =  (yj  +  y^) 

SIGZ 

=  1.0 

SIGZ  =(y2  +  y4) 

SIG3 

=  1.0 

SIG3  =  (yj  +  y^) 

74)^ 


10.  Locate  point  5  and  calculate  and 


L^  =  4-  SIGl  j  [vpj  +  Vpj^']  (xj  -  X4)  +  [Upj  +  Upj^l  (y^  - 

h  =  -g-  {  ['^pZ  +  ^p24\  ^^4  -  *2>  L'^pZ  +  '^pZ4^  <yz  “ 
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Pp4«'  =  [ 

+  p4»)  = 


~  ^pt  ~  ^2 

(Lj  +  L2) 

♦pZ-  t'>p2  +  ''p4] 


] 


*  =  r  *  P'*  ' 

L  n.  -v 

75  =  72  t  ^  (7,  -  y2> 

X5  =  +  t  (Xj  -  X^) 

By  use  of  the  above  steps,  point  5  is  located  and  and  ^  are  calculated. 
This  procedure  must  be  repeated  for  each  particle. 


11.  Evaluate  the  gas  properties  at  point  4: 

u=  ^  (Pg3  *  Pg34>  <"g3  +  "g34> 

V  =  4  (Pg3  t  Pg34)  l''g3  *  ''g34> 

''i '  T  [  'V  *  "e‘4’  ‘’'4  -  '^i*  ■  •''gi  *  ''gn’  ‘*4  -  *1>] 

‘‘ll  =  K  '“g2  *  "gZ^*  •’'4  ■  ’'2>  -  •''gZ  *  ''g24>  <*4  -  *2*  ] 

Mj  =  [  16.069  YgR  (Tg,  +  Tg,g)  (Vg  -  y,)  -  ^  K,  (u^i  4  ] 

Mjj  =  [  16. 069  YjR  (Tj2  +  T^24'  '’'4  ’  ''g’  ’  T  ‘'u  '"g2  +  "g24l] 
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Qj  4.  0174  Y  gR  (Tgi  +  (x^  -  x^) 

Qjj  =  4.  0174  VgR  (Tg2  +  gZ^)  (Pgz  +  PgZ4^  t'^gZ  ■*■  ''g24^  ^*4  "  ^2^ 

Rl  =-4.0174  YgR  (Tg^  +  (pg^  +  Pgi4'(UgjfUgj4)(x4  -  x^) 

Rjl  =-4.  0174  YgR  (Tg2  +  Tg24)  (pg2  +  Pg24)^"g2  ^  '^g24)  ^*4  "  *2^ 

6 

D  =  -  :?  I  {  A  {Pp3  +  Pp34)  [  (Ug3  +Ug34  .Up3  (x^  -  X3) 

i=l 

*  <''g3  +  ''g34  -  ''p3  -  •''4  -  ''3>  ]  } 

®I  '  [  <“gl  ♦  "gl4  -  “pi  -“pm'^  +  *''gl^''gl4  -  “pi  •  “pu)^ 

4  1.33333  C  (Tp,  4Tp,,-Tg,  -Tg,^)] 

Wg*  ^ )  p  /  2 

®n  =  -t -  L  ‘“g2  *  "g24  -  "p2  -  "p24>  +  (“gZ  *  '’g24  '  V  '  V't> 

4l.33333  C(Tp2  4Tp,,.T^2-V4>] 

(Y,-t)  r  2 

®C  “  -I -  L  '"g3  *  "g34  -  “p3  -  “p34>  *  '“g3  +  ’'g34  '  “pZ  '  ''p}4> 

+  1.33333  C(Tp3  4  Tp3^.T^,.T^3^)  ] 


2 
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h-  <'>gl"'>gl4"V  +  V4'^g'''V'8l4)Kl 


6 

7  Kj  [  [  A  Bj(pp,  4  Pp,^)  ]  .4.  0174  VgR(Tgi+ 

i=l 


I  (Ppl  +  Ppl4)  {(^gl  +  ^gl4  -  ^pl  -  ^pl4)  <^4  ■ 

i=l 


-  <''gl  +  'gU  -  ''  pi  -  ''pl4>  <*4 


-*l'}| 


(^4  -  Xj) 


=  {l^gKSlO  ,p^^  4  p^2,)  (v^,  4  VgR  (Tg,  4  T^24)  Kq 

6 

-7  KlI  Z  [a  Bjilp  24p  )  ]  -4.0174V  R(T  4 


6 

I  A  <>‘p2  *  I‘p24>  {  '"g2  “g24  '  V  '  '•p24>  '  ^2* 

i=l 


•  <''g2  *  ''g24  -  V  -  Va’  ‘*4  ■  *2>} 


<*4  •  ’‘gl 


y  =  Sj  4  4,  628  M,  (Pg,  -  Pgj)  4  Q,Ug,  4  RjV^j 


Z  =  S^  4  4,  628  Mjj  (Pg2  -  Pgj)  4  QjjUg2  +  "^n''g2 


X.D40u^,4Vv^, 


-  2 


'  a3  ^  1=1 


fg'*  g3"^g34'  '“g3  ^  “g34' 

DENOM  =[{Qj-MjTr)  (Rjj  -  MjjTT)  -  (Q^j  .  M^7)  (Rj  -  M^V)  ] 
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“,4  =  [ 


(Y  -  MjX)  (Rjj  -  MjjV)  -  (Z  -  MjjX)  (Rj  -  MjV) 


DENOM 


^g4  “  [ 


(Qj  -MjIT)  (Z  -  MjjX)  -  (Qjj  -  MjjU)  (Y  -  MjX) 


DENOM 


^g4  =  Pg3  +  2-  160  X  10  (X  -  rrug4  -  VVg^) 


288.0  (P  P) 

Pg4"  Pg3'^'^-^  [  YgR(Tg3  +  Tg34)  ] 


144  P 


12.  Evaluate  the  particle  properties  at  point  4. 


(vi  +  u  c  "  u  >«c  “  u  e) 

X=A  L  ^45  + 


t(v  ^  +  V  c  -  V  .c  -  V  c) 

g4  g5  p45  p5^  1 

p45  pS^  “* 


0-  2.6667  x  10-5  AC  [  ^  V  '  "^^4  '  j 

^  111  X  11  \  J 


<  Vs  ^  v^ 


V  "  '"pS  +  ^  <*4  ■  *5) 


V  "  V"  ®  ‘^4  -  ^5) 


V  =  V  ^  ^  ^"^4  -  *5) 


Tp4  =  ^  ( V> 
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13.  Thus,  the  location  of  point  4,  and  all  the  flow  properties  at  point  4, 
are  known  approximately.  To  improve  the  accuracy  of  the  approximation, 
point  4  must  be  relocated,  based  on  the  average  values  of  flow  properties,  by 
using  the  last-determined  values  of  the  flow  properties  at  point  4,  Then,  all 
the  flow  properties  can  be  recalculated  by  using  these  same  average  values. 

14.  Repeat  step  6  to  relocate  point  4.  Repeat  step  7  to  relocate  point  3. 

Repeat  step  8  to  set  all  as  just  determined.  Repeat  step  9  to  recal¬ 

culate  SIG  1,  SIG  2,  and  SIG  3.  Repeat  step  10  to  recalculate  Pp^j  ^  p4 

and  to  locate  point  5.  Repeat  step  11  to  recalculate  gas  properties  at  point  4. 
Repeat  step  12  to  recalculate  particle  properties  at  point  4.  Thus,  point  4  is 
relocated,  and  all  the  flow  properties  at  point  4  have  been  recalculated.  Apply 
the  following  convergence  test  on  all  flow  properties: 


P  (n  +  1)  -  Pj’n) 

- PUr+TY 


^  L. 


When  the  test  has  been  satisfied  by  all  flow  properties,  the  solution  for  the  loca¬ 
tion  and  properties  at  point  4  is  complete.  If  the  solution  does  not  converge,  the 
entire  procedure  must  be  repeated  until  convergence  is  achieved. 
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EXAMPLE  CALCULATION  FOR  VERIFI  ZCATION  OF  COMPUTER  PROGRAM 


To  check  the  general  validity  of  botbo  the  theory  and  accuracy  of  the  computer 
program,  a  numerical  example  was  calcuMlated  for  a  general  interior  point.  A 
solid  propellant  that  contains  aluminum  w^5$  chosen.  The  exhaust  gases  contained 
0.  48  lb  of  aluminum  oxide  per  lb  of  gas. 

Particle  properties  were  defined  lisx-ting  six  discrete  particle  sizes  (Table  I- 1). 

The  transport  properties  of  the  gas  phase  atchsLxnber  stagnation  conditions  were 

defined  and  are  given  in  Table  1-2.  T.  v _ ,  and  C  are  given  as  functions  of  gas 

8  PS 

temperature  in  Table  1-3.  The  temperatuzcre-enthalpy  relationship  for  solid  and 
liquid  aluminum  oxide  is  shown  in  Table  I — 4  with  other  particle  properties  needed 
for  the  solution.  The  f  and  g  parameters  a  .8 functions  of  Reynolds  number  were  ob¬ 
tained  from  current  studies  and  are  presetmted  in.  Table  1-5.  Values  listed  in  this 
table  may  be  considered  ap  permanent  valujues  of  f  and  g  for  all  subsequent  work. 
Units  of  physical  parameters  are  shown  in  Table  1-6. 

The  general  and  specific  data,  outliaeedin  .Appendix  H,  are  shown  in  Tables 
1-7  and  -8,  respecitvely.  The  example  cal»alation  is  based  on  the  characteristic  net 
shown  in  Figure  I-l  and  the  location  and  pcaroperti es  for  the  general  interior  point  4 
are  given  in  Table  1-9. 

The  gas  viscosity  for  several  propelMlants  was  found  to  satisfy  the  relationship. 


Hence,  if  \l^  is  known  at  any  temperature, 


fi  may  be  determined.  Usually  jX  is 
®  g 


known  at  the  flame  temperature  as  a  results  of  the  rmochemical  calculations. 


In  comparison  of  results,  u  and  v  cr-decreased,  and  p  and  P  remained 

8  8  8  8 

about  the  same,  and  T^  increased,  which  iss  almost  opposite  to  the  trend  expected 
in  a  perfect-gas  flow.  However,  in  the  exa::.~inple  calculation, u^  and  v^  increased, 
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and  pp  and  decreased,  which  was  the  expected  trend.  The  apparent  inconsist¬ 
encies  in  the  gas  properties  are  due  to  the  fact  that  this  problem  was  artificially 
conceived.  An  actual  problem  beginning  at  the  supersonic  starting  line  should 
behave  in  a  more  normal  manner. 


The  example  calculation  appears  to  verify  the  theory  developed  in  the  pre- 
ceeding  appendixes  and  authenticate  the  correctness  of  the  computer  program. 
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Table  I-l 


PARTICLE  PROPERTIES 
Chemical  Species  =  AI2O3 


Mass  flow  rate  of  particles 
Mass  flow  rate  of  gas 

Particle  Radii,  microns 


"■pl  " 
'■p3  " 

^p4  " 
^p6  = 


0.45 
0.  72 
0.90 
1.  15 
1.42 
1.75 
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Table  1-2 

GAS  TRANSPORT  PROPERTIES 


Pressure,  P 

g 

1 , 000. 0  psia 

Temperature,  Tg 

6,650.0  ‘R 

Molecular  Weight  (Gas) 

20.99 

Gas  Constant,  R 

o  ft-lbf 

Ibrn-'R - 

Specific  Heat  Ratio,  Y  g 

1.20 

Constant  Pressure  Specific  Heat,  c 

Pg 

0.4776 

Btu 

'  Ibrn-’R  ■ 

Viscosity,  p, 

g 

0.2416 

Ibm 

Density,  p 

0.2940 

Ibm 

g 

ft3 

Prandtl  Number,  Pr 

0.4996 

Thermal  Conductivity, 

0.2307 

Btu 

Tt-hr-^R 
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200 
500 
1,  000 
2,  000 

3,  000 

4,  000 

5,  000 

6,  000 

6.  650 


Table  1-3 


AND  C  VS  TEMPERATURE 
Pg 


Btu 

'lbm:*R 


1.366 

1.365 

1.344 

1.297 

1.257 

1.244 

1.228 

1.212 

1.200 


0.2860 
0.2870 
0. 3130 
0.3730 
0.4210 
0.4520 
0.4700 
0.4750 
0.4776 
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Table  1-4 

PROPERTIES  OF  ALUMINUM  OXIDE 


T  =4,  172.4  ‘R 
pm 


1,  051.4 


Btu 

Ibm 


1,695.5 


Btu 

Ibm 


Solid 


Liquid 


,  Btu 

p’  Ibm 

Tp.  -R 

-42.26 

0 

-40.89 

180 

0.62 

540 

37.90 

720 

80.80 

900 

175.5 

1,260 

276.9 

1,620 

435.6 

2,  160 

598.6 

2,700 

874.7 

3,600 

1,  051.4 

4,  172.4 

Btu 
p’  Ibm 

Tp.  «R 

1,695.5 

4,  172.4 

1,745.6 

4,320 

1,8Q7.2 

4,500 

3,  043.0 

8,  100 
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Table  1-5 

f  AND  g  PARAMETERS  AS  FUNCTIONS  OF  REYNOLDS  NUMBER  * 


Re 

f 

Re 

g 

0 

1.0000 

0 

1. 0000 

1 

1,0000 

3 

1. 0000 

4 

1.4167 

10 

1. 3750 

10 

1.7917 

20 

1.7000 

40 

3. 1667 

30 

2, 0000 

70 

4. 0800 

40 

2. 1500 

100 

4.7917 

70 

2.8000 

200 

6.6667 

100 

2.9300 

500 

12. 1000 

1,  000 

16.6670 

For  Re  ^  1, 000, 

For  Re  ^  100, 

f  =  0.016667  Re 

g  =  0.  185  (Re) 

*  f  data  taken  from  Boundary  Layer  Theory  by  Schiichting 
g  data  taken  from  Heat  Transmission  by  McAdams 
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TABLE  1-6 

UNITS  OF  PHYSICAL  PARAMETERS 


specific  heat  at  constant  pressure 
specific  heat  ratio  of  the  gas 
ratio  of  Nu  to  Nu  for  Stokes'  flow 
ratio  of  Nu  to  Nu  for  Stokes'  flow 
particle  enthalpy 

particle  density  per  unit  volume  of  particle 
gas  static  pressure 
Prandtl  number 

particle  radius,  baaed  on  spherical  particles 

gas  constant 

Reynolds  number 

gas  static  temperature 

particle  static  temperature 

axial  gas  velocity  component 

axial  particle  velocity  component 

normal  or  radial  gas  velocity  component 

normal  or  radial  particle  velocity  component 

axial  coordinate 

component  of  gas  viscosity  expression 
gas  viscosity 

coefficient  of  gas  viscosity  expression 
gas  density 

particle  density  per  unit  volume  of  gas 


Btu/lbm-  ®R 

dimensionless 

dimensionless 

dimensionless 

Btu/lbm 

Ibm/  cu  ft 

Ibf/sq  in. 

dimensionless 

micron 

Btu/lbm-  •R 

dimensionless 

•R 

•R 

ft/aec 

ft/sec 

ft/sec 

ft/sec 

ft 

dimensionless 
Ibm/ft-sec 
lbm/ft-sec(  ‘R)** 
Ibm/cu  ft 
Ibm/cu  ft 
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TABLE  1-7 

GENERAL  DATA 


Axisymmetric  flow,  =  1 

Six  aluminiom -oxide  particles,  m  s  247. 61,  with  the  radii  given  in 

Table  I- 1.  P 

f  is  given  in  Table  1-5 

g  is  given  in  Table  1-5 

Y  is  given  in  Table  1-3 

c  is  given  in  Table  1-3 
Pg 

T  vs  h  is  given  in  Table  1-4 
P  P 

T  =4,  172.4 

pm 

Pr  =  0.4996 

u  =  6.6  x  10~®  and  a  -  0.787 
o 


R  =  72.8 
K  =  L  =  0.  0001 
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TABLE  1-8 


SPECIFIC  DATA  AT  INITIAL  POINTS  1  AND  2* 

1.  =  1.00000 

=  0.60000 
2-  =  8,960.0 


Vgi  =  2.460.0 
Pgj  =  0. 00022444 
Pgl  =  0.45500 
Tgj  =  4,  010.  0 
'^pl  (1)  =  8,  050 
Vpi  (1)  =  2,  140 
Ppj  (1)  =  0. 0000195 
hpi  (1)  =  1. 035.2 
Tpi(l)  =  4,120.0 


"pi 

(2)  : 

=  7,  650 

'"pi 

(2)  = 

=  1, 880 

•^pl 

(2) 

=  0. 0000862 

h  , 
pi 

(2)  : 

=  1, 054. 0 

T  , 
pi 

(2) 

=  4,  172.4 

"pi 

(3)  = 

=  7,220 

'"pl 

(3)  = 

1,640 

Ppl 

(3)  = 

0. 0000502 

^1 

(3)  = 

1, 300. 0 

^Pl 

(3)  = 

4.  172.4 

^2  =  1. 02500 
72  -  0.50000 
Ug2  =  9,  060.  0 
Vg2  =  2.  370.0 
Pg2  =  0. 00018289 
Pg2  0.  36800 
Tg2  =  3,  980.  0 

Up2  (1)  ::  8,  100 
Vp2  (1)  =  2.200 

Pp2  (1)  -  0. 0000160 
hp2  (1)  =  1, 029. 1 
Tp2  (1)  =  4,  100.0 

V  =  7.710 

V  f2)  =  1. 940 
Pp2  (2)  =  0.  0000710 

V  052.0 

^p2  172.4 

'^p2  =  7.280 

^p2  <3)  =  1.  700 

Pp2  =  0-  0000412 

V  =  1.295.0 

'^p2  =  4.  172.4 
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TABLE  1-8  (cont.) 


(4)  =  6.790 
(4)  =  1,420 
Ppl  (4)  =  0.  0000302 
hp^  (4)  =  1, 696. 0 
Tp^  (4)  =  4,  173.9 

Up^  (5)  =  6.  380 
Vp^  (5)  =  1.225 
Pp^  (5)  =  0. 0000211 
hp^  (5)  =  1,697.3 
Tpj  (5)  =  4,  177.7 

(6)  =  5.930 
Vpl  (6)  =  1,030 
Ppl  (6)  =  0. 0000118 
hpi  (6)  =  1,779.8 
Tpi  (6)  =  4,  420.  0 


Up2  (4)  =  6,850 
Vp2  (4)  =  1, 480 
Pp2  (4)  =  0.  0000249 
hp2  (4)  =  1, 694.  5 
Tp2  (4)  =  4,  172.4 

Up2  (5)  =  6,  420 
Vp2(5)=  1.265 
Pp2  (5)  =  0. 0000173 
hp2  (5)  =  1, 696.  0 
Tp2  (5)  =  4,  173.9 

V  =  5.970 

V  (6)  =  1,  070 

Pp2  ^5)  =  0.  OOOOr'O? 
^p2  773. 0 

Tp2  (6)  =  4,400.  0 


*  The  location  and  properties  of  point  4  (Figure  I-l  and  Table  1-8)  for  the  above 
case  are  given  in  Table  1-9  in  the  format  printed  out  by  the  computer. 
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